Acknowledgments

During the preparation of this work, we had discussions with the group formed by Elie Ghazzoul and Matteo Toschi to address doubts encountered while solving the problems. Additionally, we declare that ChatGPT was used as support during debugging and for verifying the code.

Part I

Point a

The following section provides a formal description of conformal prediction and adaptive conformal inference (ACI). The content is based on the paper by Gibbs and Candes (2021).


Conformal Inference Basics (Section 1.1)

Conformal prediction is a technique in machine learning that transforms any prediction model into one that outputs prediction intervals or sets. These intervals are unique because they come with a formal guarantee: the true value of the target variable will lie within the interval a specific percentage of the time, such as 90% or 95%.

This is revolutionary because:

  1. It works with any machine learning model, whether it’s a simple linear regression or a complex neural network.

  2. It gives reliable and understandable results, which is critical for high-stakes applications like finance or healthcare.

  3. It quantifies uncertainty, making black-box models more transparent and trustworthy.

The ability to provide guaranteed coverage is what sets conformal prediction apart from traditional methods, which often struggle to deliver reliable uncertainty estimates, especially for unseen data.

How Does Conformal Prediction Work?

Conformal prediction relies on three key components: conformity scores, the quantile function, and the prediction set. Let’s break these down.


1. Conformity Score

The conformity score is the heart of the conformal prediction framework. It measures how well a candidate prediction aligns with the model’s output. In simple terms, it tells us if a prediction is “reasonable” or not, based on past data.

  • For a regression model predicting a single value \(\hat{\mu}(X)\) for a given input \(X\), a natural conformity score is:

    \[ S(X, Y) = |\hat{\mu}(X) - Y| \]

    This measures the absolute difference between the predicted value \(\hat{\mu}(X)\) and the actual observed value \(Y\). Smaller scores mean the prediction is closer to the true value.

  • For models that predict quantiles instead of single values (e.g., upper and lower bounds), a more sophisticated conformity score might be:

    \[ S(X, Y) = \max\{\hat{q}(X; \alpha/2) - Y, Y - \hat{q}(X; 1 - \alpha/2)\} \]

    Here, \(\hat{q}(X; p)\) represents the predicted quantile for a given probability \(p\). This score checks how far \(Y\) is from the predicted upper or lower quantiles.

Why is this important?
The conformity score allows us to compare the new prediction against the model’s past performance, making it possible to decide whether the prediction is acceptable.


2. Quantile Function

The quantile function is how we decide the threshold for what counts as a “reasonable” prediction. It tells us how large the conformity score can be while still keeping the prediction within a desired level of certainty.

To calculate the quantile function, we use a calibration dataset \(D_\text{cal}\), which is separate from the data used to train the model. The calibration dataset contains past data points \((X, Y)\), and we compute the conformity scores for these points. Then, we find the quantile \(Q(p)\) such that:

\[ Q(p) = \inf \left\{s : \frac{1}{|D_\text{cal}|} \sum_{(X, Y) \in D_\text{cal}} \mathbf{1}[S(X, Y) \leq s] \geq p \right\} \]

In simple terms:

  • The quantile \(Q(p)\) is the smallest value of the conformity score \(S\) that ensures at least \(p\) fraction of the calibration scores are below it.

  • For example, if \(p = 0.9\), \(Q(0.9)\) will capture the 90th percentile of conformity scores in the calibration dataset.

Why is this important?
The quantile function dynamically adjusts the prediction intervals based on past data. This ensures that the intervals remain valid and reflect the actual uncertainty in the data.


3. Prediction Set

Finally, we use the conformity score and quantile function to construct the prediction set, which is the interval we provide for the new prediction.

For a new input \(X_t\), the prediction set \(\hat{C}_t\) is defined as: \[ \hat{C}_t = \{y : S(X_t, y) \leq Q(1 - \alpha)\} \]

  • Here, \(\alpha\) is the desired error rate (e.g., \(\alpha = 0.1\) for 90% coverage).

  • The set contains all possible values \(y\) that have a conformity score \(S(X_t, y)\) below the threshold \(Q(1 - \alpha)\).

Why is this important?
The prediction set provides a range of values that are most likely to include the true value \(Y_t\). This range is guaranteed to achieve the desired coverage level, assuming the data meets the exchangeability condition.


Adaptive Conformal Inference (Section 2)

Traditional conformal prediction methods assume that the data is *exchangeable, meaning that the order of the data points does not matter, and the distribution of the data remains constant over time. However, this assumption often fails in real-world applications, especially for **time series data* or other non-stationary data where the data distribution can shift significantly over time.

Adaptive Conformal Inference (ACI) addresses this challenge by dynamically recalibrating the coverage parameter \(\alpha_t\) as the data evolves. This makes ACI robust to changes in the data’s distribution, ensuring that the prediction intervals maintain their validity over time.


Key Concepts of Adaptive Conformal Inference

1. Miscoverage Rate

The miscoverage rate quantifies how often the true value \(Y_t\) falls outside the prediction set \(\hat{C}_t\).

For a given time step \(t\), the miscoverage rate \(M_t(\alpha)\) is defined as: \[ M_t(\alpha) = P(S(X_t, Y_t) > Q_t(1 - \alpha)) \] This expression states that the miscoverage rate is the probability that the conformity score \(S(X_t, Y_t)\) exceeds the threshold \(Q_t(1 - \alpha)\), where \(Q_t(1 - \alpha)\) is the quantile of the conformity scores determined from past data.

In practical terms:

  • If \(M_t(\alpha)\) is close to \(\alpha\), the prediction intervals are well-calibrated.

  • If \(M_t(\alpha)\) is far from \(\alpha\), the intervals are either too narrow (under-covering) or too wide (over-covering).


2. Adaptive Update Rule

The core innovation of ACI is its adaptive update rule, which recalibrates the coverage parameter \(\alpha_t\) at each time step based on observed performance.

The update rule is: \[ \alpha_{t+1} = \alpha_t + \gamma (\alpha - \text{err}_t) \]

Explanation of the Terms:

  1. Current Parameter (\(\alpha_t\)): The coverage parameter at time \(t\).

  2. Step Size (\(\gamma\)): A small, positive constant that controls how quickly the method adapts to changes in the data. Larger \(\gamma\) allows for faster adaptation but may introduce volatility, while smaller \(\gamma\) provides stability at the cost of slower adaptation.

  3. Target Coverage (\(\alpha\)): The desired miscoverage rate (e.g., \(\alpha = 0.1\) for 90% coverage).

  4. Error Term (\(\text{err}_t\)):

    • \(\text{err}_t = 1\) if the true value \(Y_t\) is outside the prediction set \(\hat{C}_t\).

    • \(\text{err}_t = 0\) if the true value \(Y_t\) is inside the prediction set \(\hat{C}_t\).

This update rule works by:

  • Increasing \(\alpha_t\) if the prediction set has been under-covering (too narrow, too many \(\text{err}_t = 1\)).

  • Decreasing \(\alpha_t\) if the prediction set has been over-covering (too wide, too few \(\text{err}_t = 1\)).


3. Objective of ACI

The primary objective of ACI is to ensure that the coverage level \(1 - \alpha\) is maintained approximately over time, even as the data distribution shifts. By dynamically adjusting \(\alpha_t\), ACI can respond to non-stationarity in the data while maintaining reliable prediction intervals.


ACI in Practice

Here’s how ACI works step by step: 1. Compute the Conformity Score:

  • For each new data point \(X_t\), calculate the conformity score \(S(X_t, Y_t)\) for a candidate prediction \(Y_t\).
  1. Construct the Prediction Set:

    • Use the current \(\alpha_t\) to determine the prediction set: \[ \hat{C}_t = \{ y : S(X_t, y) \leq Q_t(1 - \alpha_t) \} \]

    • The set \(\hat{C}_t\) contains all values \(y\) that are considered reasonable predictions based on the conformity score and the threshold \(Q_t(1 - \alpha_t)\).

  2. Evaluate Miscoverage:

    • Check whether the true value \(Y_t\) falls within the prediction set \(\hat{C}_t\):

      • If \(Y_t \notin \hat{C}_t\), record a miscoverage event (\(\text{err}_t = 1\)).

      • Otherwise, record \(\text{err}_t = 0\).

  3. Update \(\alpha_t\):

    • Apply the adaptive update rule to adjust \(\alpha_t\) based on recent performance.

Why Is ACI Important?

  1. Handles Non-Stationary Data:

    • Unlike traditional conformal methods, ACI is designed for data where the distribution can change over time, such as financial markets or environmental data.
  2. Simple and Efficient:

    • The update rule is computationally lightweight and only requires tracking a single parameter (\(\alpha_t\)) over time.
  3. Robust Prediction Intervals:

    • ACI ensures that the prediction intervals remain reliable even under significant distribution shifts, providing greater trust in the predictions.

Point b: Implementation and Comparison of Strategies

In this section, we applied and compared different conformal prediction strategies to quantify the prediction uncertainty of Apple stock volatility. Following the instructions in the assignment, we implemented the Standard Conformal Prediction (CP), the Adaptive Simple strategy (based on Equation 2), and the Momentum Adaptive strategy (based on Equation 3) using both GARCH(1,1) and Exponential GARCH (EGARCH)(1,1) models. The main difference between EGARCH and GARCH models is that EGARCH accounts for asymmetries in volatility, allowing negative shocks to have a different impact on volatility compared to positive shocks, while GARCH assumes a symmetric response to all shocks.

Methodology

  1. Target Parameters:

    • The desired coverage level was set to \(90\%\) (\(\alpha = 0.1\)).

    • For the adaptive strategies:

      • The Adaptive Simple strategy and Momentum Adaptive strategy were initially tested with an adaptation rate \(\gamma = 0.05\).
  2. Conformity Scores:

    • For Standard CP and Adaptive Simple, we used both normalized and non-normalized conformity scores.

    • For Momentum Adaptive, only normalized conformity scores were used.

  3. Data:

    • Apple stock data from 2015-01-01 to 2023-12-31 was used.

    • Both GARCH and Exponential GARCH models were fitted to the data to generate volatility forecasts.

  4. Evaluation of Strategies:

    • Prediction intervals were constructed for each method, and local coverage rates were calculated to assess whether the observed volatility fell within the predicted intervals at the target \(90\%\) level.

    • Additional plots were created to track the behavior of prediction errors over time, ensuring they converged to the expected \(10\%\) miscoverage rate (\(\alpha = 0.1\)).

  5. Model Selection with \(\gamma = 0.05\):

    • For each strategy and model, the local coverage rates were compared against the target level.

    • The best-performing model was selected by minimizing the sum of the absolute difference between the mean local coverage rate and the target, and the variance of the local coverage rates.

  6. Testing Different Values of \(\gamma\):

    • To evaluate the impact of the adaptation rate \(\gamma\), we repeated the analysis using \(\gamma = 0.03\) and \(\gamma = 0.07\), focusing only on the three best-performing strategies identified with \(\gamma = 0.05\).

    • For each \(\gamma\), the local coverage rates and prediction errors were plotted, and the best-performing strategy was selected.

    • Finally, the best strategies for \(\gamma = 0.03\), \(\gamma = 0.05\), and \(\gamma = 0.07\) were compared to determine the overall best-performing model and adaptation rate.

Objective

The main objective of this analysis was to compare the coverage performance of the different strategies and evaluate how the adaptation rate \(\gamma\) affects the results. The final goal was to identify the strategy and model that consistently provide reliable prediction intervals at the target coverage level.


# Load necessary libraries
library(quantmod)
library(quantreg)
library(rugarch)
library(dplyr)
library(zoo)

### Volatility prediction
#Garch
garchConformalForcasting <- function(returns, alpha = 0.1, gamma = 0.05, lookback = 1250, 
                                     garchP = 1, garchQ = 1, startUp = 100, verbose = FALSE, 
                                     updateMethod, momentumBW = 0.8, scoreType ) {
  myT <- length(returns)
  T0 <- max(startUp, lookback)
  garchSpec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), 
                          variance.model = list(model = "sGARCH", garchOrder = c(garchP, garchQ)), 
                          distribution.model = "norm")
  alphat <- alpha
  ### Initialize data storage variables
  errSeqOC <- rep(0, myT - T0 + 1)
  errSeqNC <- rep(0, myT - T0 + 1)
  alphaSequence <- rep(alpha, myT - T0 + 1)
  scores <- rep(0, myT - T0 + 1)
  
  for (t in T0:myT) {
    if (verbose) {
      print(t)
    }
    ### Fit GARCH model and compute new conformity score
    garchFit <- ugarchfit(garchSpec, returns[(t - lookback + 1):(t - 1)], solver = "hybrid")
    sigmaNext <- sigma(ugarchforecast(garchFit, n.ahead = 1))
    if (scoreType == "Simple") {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2)
    } else {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2) / sigmaNext^2
    }
    
    recentScores <- scores[max(t - T0 + 1 - lookback + 1, 1):(t - T0)]
    
    ### Compute errt for both methods
    errSeqOC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alphat))))
    errSeqNC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alpha))))
    
    ### Update alphat
    alphaSequence[t - T0 + 1] <- alphat
    if (updateMethod == "Simple") {
      alphat <- alphat + gamma * (alpha - errSeqOC[t - T0 + 1])
    } else if (updateMethod == "Momentum") {
      w <- rev(momentumBW^(1:(t - T0 + 1)))
      w <- w / sum(w)
      alphat <- alphat + gamma * (alpha - sum(errSeqOC[1:(t - T0 + 1)] * w))
    }
    if (t %% 100 == 0) {
      print(sprintf("Done %g steps", t))
    }
  }
  
  return(list(alphaSequence = alphaSequence, errSeqOC = errSeqOC, errSeqNC = errSeqNC))
}

# Step 1: Get financial data
getSymbols("AAPL", from = "2015-01-01", to = "2023-12-31")
## [1] "AAPL"
aapl_prices <- Cl(AAPL)
returns <- dailyReturn(aapl_prices)
returns <- na.omit(returns)

results_Studentized <- garchConformalForcasting(returns = returns, alpha = 0.1, gamma = 0.05,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
results_Simple <- garchConformalForcasting(returns = returns, alpha = 0.1, gamma = 0.05,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Simple")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
date <- index(aapl_prices)[1250:length(index(aapl_prices))]
alphaSequence_St <- results_Studentized[[1]]
errSeqOC_St <- results_Studentized[[2]]; mean(errSeqOC_St)
## [1] 0.09852217
errSeqNC_St <- results_Studentized[[3]]; mean(errSeqNC_St)
## [1] 0.08078818
alphaSequence_Si <- results_Simple[[1]]
errSeqOC_Si <- results_Simple[[2]]; mean(errSeqOC_Si)
## [1] 0.09162562
errSeqNC_Si <- results_Simple[[3]]; mean(errSeqNC_Si)
## [1] 0.06699507
results_StudentizedM <- garchConformalForcasting(returns = returns, alpha = 0.1, gamma = 0.05,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Momentum",scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
alphaSequence_StM <- results_StudentizedM[[1]]
errSeqOC_StM <- results_StudentizedM[[2]]; mean(errSeqOC_StM)
## [1] 0.09753695
errSeqNC_StM <- results_StudentizedM[[3]]; mean(errSeqNC_StM)
## [1] 0.08078818
plot(date, cummean(errSeqNC_St), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: GARCH Strategies Studentized", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date, cummean(errSeqOC_St), col = "blue")
lines(date, cummean(errSeqOC_StM), col = "red",type = "s")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple", "Adaptive Momentum"), 
       col = c("black", "blue", "red"), lty = 1)

plot(date, cummean(errSeqNC_Si), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: GARCH Strategies Simple", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date, cummean(errSeqOC_Si), col = "blue")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple"), 
       col = c("black", "blue"), lty = 1)

compute_local_coverage <- function(errors, window = 250) {
  1 - rollmean(errors, k = window, fill = NA)
}

window_size <- 250

# Generate Local Coverage Rates for Normalized and Non-Normalized Scores

# Local Coverage for Scores Non-Normalized
local_coverage_fixed_non_normalized <- compute_local_coverage(errSeqNC_Si, window = window_size)
local_coverage_adaptive_simple_non_normalized <- compute_local_coverage(errSeqOC_Si, window = window_size)

# Local Coverage for Scores Normalized
local_coverage_fixed_normalized <- compute_local_coverage(errSeqNC_St, window = window_size)
local_coverage_adaptive_simple_normalized <- compute_local_coverage(errSeqOC_St, window = window_size)
local_coverage_adaptive_momentum <- compute_local_coverage(errSeqOC_StM, window = window_size)

# Plot Local Coverage Rates for Non-Normalized Scores
plot(date, local_coverage_adaptive_simple_non_normalized, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "GARCH Local Coverage Rates: Non-Normalized", xlab = "Date", ylab = "Local Coverage Level")
lines(date, local_coverage_fixed_non_normalized, col = "red")
abline(h = 0.9, col = "black", lty = 2)
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha"), 
       col = c("blue", "red"), lty = 1)

# Plot Local Coverage Rates for Normalized Scores
plot(date, local_coverage_adaptive_simple_normalized, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "GARCH Local Coverage Rates: Normalized", xlab = "Date", ylab = "Local Coverage Level")
lines(date, local_coverage_fixed_normalized, col = "red")
lines(date, local_coverage_adaptive_momentum, col = "green", type="s")
abline(h = 0.9, col = "black", lty = 2)  # Target coverage level
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha", "Adaptive Momentum"), 
       col = c("blue", "red", "green"), lty = 1)

Results: Error Rates (Studentized)

For the studentized conformity scores, we observe that both the Adaptive Simple and Adaptive Momentum strategies closely approach the target error rate (\(\alpha = 0.1\)), while the Standard CP shows a greater distance from the target. Additionally, we can see that the performance of Adaptive Simple and Adaptive Momentum is nearly identical, with both strategies achieving very similar results.


Results: Error Rates (Non-Normalized)

When using non-normalized conformity scores, there is a noticeable difference between the Adaptive Simple strategy and Standard CP. The Adaptive Simple strategy demonstrates a much better convergence toward the target error rate, while the Standard CP remains further from the target throughout the evaluation period.


Results: Local Coverage Rates (Non-Normalized)

For non-normalized conformity scores, we observe a clear graphical difference between the Adaptive Simple and Standard CP strategies. While both strategies oscillate around the target coverage level (\(90\%\)), the Standard CP tends to remain significantly higher than the target, whereas the Adaptive Simple strategy consistently converges more closely to the target.


Results: Local Coverage Rates (Normalized)

For normalized conformity scores, the difference between strategies becomes even more pronounced. The Standard CP remains above the target coverage level for almost the entire evaluation period, failing to reach the desired level. In contrast, both the Adaptive Simple and Adaptive Momentum strategies achieve excellent performance, closely fluctuating around the target coverage level and showing very similar results overall.

### eGARCH
EgarchConformalForcasting <- function(returns, alpha = 0.1, gamma = 0.05, lookback = 1250, 
                                      garchP = 1, garchQ = 1, startUp = 100, verbose = FALSE, 
                                      updateMethod, momentumBW = 0.8, scoreType ) {
  myT <- length(returns)
  T0 <- max(startUp, lookback)
  
  # Updated to eGARCH specification
  garchSpec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), 
                          variance.model = list(model = "eGARCH", garchOrder = c(garchP, garchQ)), 
                          distribution.model = "norm")
  
  alphat <- alpha
  errSeqOC <- rep(0, myT - T0 + 1)
  errSeqNC <- rep(0, myT - T0 + 1)
  alphaSequence <- rep(alpha, myT - T0 + 1)
  scores <- rep(0, myT - T0 + 1)
  
  for (t in T0:myT) {
    if (verbose) {
      print(t)
    }
    # Fit eGARCH model and compute conformity score
    garchFit <- ugarchfit(garchSpec, returns[(t - lookback + 1):(t - 1)], solver = "hybrid")
    sigmaNext <- sigma(ugarchforecast(garchFit, n.ahead = 1))
    
    if (scoreType == "Simple") {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2)
    } else {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2) / sigmaNext^2
    }
    
    recentScores <- scores[max(t - T0 + 1 - lookback + 1, 1):(t - T0)]
    errSeqOC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alphat))))
    errSeqNC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alpha))))
    
    alphaSequence[t - T0 + 1] <- alphat
    if (updateMethod == "Simple") {
      alphat <- alphat + gamma * (alpha - errSeqOC[t - T0 + 1])
    } else if (updateMethod == "Momentum") {
      w <- rev(momentumBW^(1:(t - T0 + 1)))
      w <- w / sum(w)
      alphat <- alphat + gamma * (alpha - sum(errSeqOC[1:(t - T0 + 1)] * w))
    }
    if (t %% 100 == 0) {
      print(sprintf("Done %g steps", t))
    }
  }
  
  return(list(alphaSequence = alphaSequence, errSeqOC = errSeqOC, errSeqNC = errSeqNC))
}

results_EStudentized <- EgarchConformalForcasting(returns = returns, alpha = 0.1, gamma = 0.05,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
results_ESimple <- EgarchConformalForcasting(returns = returns, alpha = 0.1, gamma = 0.05,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Simple")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
date <- index(aapl_prices)[1250:length(index(aapl_prices))]
alphaSequence_ESt <- results_EStudentized[[1]]
errSeqOC_ESt <- results_EStudentized[[2]]; mean(errSeqOC_ESt)
## [1] 0.09852217
errSeqNC_ESt <- results_EStudentized[[3]]; mean(errSeqNC_ESt)
## [1] 0.06699507
alphaSequence_ESi <- results_ESimple[[1]]
errSeqOC_ESi <- results_ESimple[[2]]; mean(errSeqOC_ESi)
## [1] 0.09261084
errSeqNC_ESi <- results_ESimple[[3]]; mean(errSeqNC_ESi)
## [1] 0.06403941
results_EStudentizedM <- EgarchConformalForcasting(returns = returns, alpha = 0.1, gamma = 0.05,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Momentum",scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
alphaSequence_EStM <- results_EStudentizedM[[1]]
errSeqOC_EStM <- results_EStudentizedM[[2]]; mean(errSeqOC_EStM)
## [1] 0.09753695
errSeqNC_EStM <- results_EStudentizedM[[3]]; mean(errSeqNC_EStM)
## [1] 0.06699507
plot(date, cummean(errSeqNC_ESt), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: eGARCH Strategies Studentized", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date, cummean(errSeqOC_ESt), col = "blue")
lines(date, cummean(errSeqOC_EStM), col = "red",type = "s")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple", "Adaptive Momentum"), 
       col = c("black", "blue", "red"), lty = 1)

plot(date, cummean(errSeqNC_ESi), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: eGARCH Strategies Simple", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date, cummean(errSeqOC_ESi), col = "blue")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple"), 
       col = c("black", "blue"), lty = 1)

compute_local_coverage <- function(errors, window = 250) {
  1 - rollmean(errors, k = window, fill = NA)
}

window_size <- 250
# Generate Local Coverage Rates for Normalized and Non-Normalized Scores

# Local Coverage for Scores Non-Normalized (eGARCH)
local_coverage_Efixed_non_normalized <- compute_local_coverage(errSeqNC_ESi, window = window_size)
local_coverage_adaptive_Esimple_non_normalized <- compute_local_coverage(errSeqOC_ESi, window = window_size)

# Local Coverage for Scores Normalized (eGARCH)
local_coverage_Efixed_normalized <- compute_local_coverage(errSeqNC_ESt, window = window_size)
local_coverage_adaptive_Esimple_normalized <- compute_local_coverage(errSeqOC_ESt, window = window_size)
local_coverage_adaptive_Emomentum <- compute_local_coverage(errSeqOC_EStM, window = window_size)

# Plot Local Coverage Rates for Non-Normalized Scores (eGARCH)
plot(date, local_coverage_adaptive_Esimple_non_normalized, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "eGARCH Local Coverage Rates: Non-Normalized", xlab = "Date", ylab = "Local Coverage Level")
lines(date, local_coverage_Efixed_non_normalized, col = "red")
abline(h = 0.9, col = "black", lty = 2)
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha"), 
       col = c("blue", "red"), lty = 1)

# Plot Local Coverage Rates for Normalized Scores (eGARCH)
plot(date, local_coverage_adaptive_Esimple_normalized, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "eGARCH Local Coverage Rates: Normalized", xlab = "Date", ylab = "Local Coverage Level")
lines(date, local_coverage_Efixed_normalized, col = "red")
lines(date, local_coverage_adaptive_Emomentum, col = "green", type="s")
abline(h = 0.9, col = "black", lty = 2)  # Target coverage level
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha", "Adaptive Momentum"), 
       col = c("blue", "red", "green"), lty = 1)

Results: Error Rates (Studentized)

For the studentized conformity scores using the eGARCH model, the results are very similar to those observed for GARCH. Both the Adaptive Simple and Adaptive Momentum strategies perform well, closely approaching the target error rate (\(\alpha = 0.1\)). The Standard CP, on the other hand, shows a greater distance from the target throughout the period. Additionally, the Adaptive Simple and Adaptive Momentum strategies once again exhibit almost identical performance.


Results: Error Rates (Non-Normalized)

For non-normalized conformity scores, the difference between Adaptive Simple and Standard CP is consistent with what was observed for GARCH. The Adaptive Simple strategy demonstrates significantly better convergence toward the target error rate, while the Standard CP remains farther from the target over time.


Results: Local Coverage Rates (Non-Normalized)

The local coverage rates for non-normalized conformity scores using the eGARCH model confirm the same pattern observed with GARCH. The Adaptive Simple strategy oscillates around the target level (\(90\%\)) but converges more closely to it compared to the Standard CP, which remains consistently higher than the target.


Results: Local Coverage Rates (Normalized)

With normalized conformity scores, the difference between strategies becomes more evident, mirroring the results of GARCH. The Standard CP strategy fails to achieve the target coverage level, remaining above it for most of the period. In contrast, both the Adaptive Simple and Adaptive Momentum strategies show excellent performance, closely fluctuating around the target and demonstrating nearly identical results.

# Compare and Select the Best Models
model_means <- list(
  # Standard CP
  Standard_CP_GARCH_NonNormalized = mean(local_coverage_fixed_non_normalized, na.rm = TRUE),
  Standard_CP_eGARCH_NonNormalized = mean(local_coverage_Efixed_non_normalized, na.rm = TRUE),
  Standard_CP_GARCH_Normalized = mean(local_coverage_fixed_normalized, na.rm = TRUE),
  Standard_CP_eGARCH_Normalized = mean(local_coverage_Efixed_normalized, na.rm = TRUE),
  
  # Adaptive Simple
  Adaptive_Simple_GARCH_NonNormalized = mean(local_coverage_adaptive_simple_non_normalized, na.rm = TRUE),
  Adaptive_Simple_eGARCH_NonNormalized = mean(local_coverage_adaptive_Esimple_non_normalized, na.rm = TRUE),
  Adaptive_Simple_GARCH_Normalized = mean(local_coverage_adaptive_simple_normalized, na.rm = TRUE),
  Adaptive_Simple_eGARCH_Normalized = mean(local_coverage_adaptive_Esimple_normalized, na.rm = TRUE),
  
  # Adaptive Momentum
  Adaptive_Momentum_GARCH = mean(local_coverage_adaptive_momentum, na.rm = TRUE),
  Adaptive_Momentum_eGARCH = mean(local_coverage_adaptive_Emomentum, na.rm = TRUE)
)

# Compare models using both mean and variance of local coverage rates
model_variances <- list(
  # Standard CP
  Standard_CP_GARCH_NonNormalized = var(local_coverage_fixed_non_normalized, na.rm = TRUE),
  Standard_CP_eGARCH_NonNormalized = var(local_coverage_Efixed_non_normalized, na.rm = TRUE),
  Standard_CP_GARCH_Normalized = var(local_coverage_fixed_normalized, na.rm = TRUE),
  Standard_CP_eGARCH_Normalized = var(local_coverage_Efixed_normalized, na.rm = TRUE),
  
  # Adaptive Simple
  Adaptive_Simple_GARCH_NonNormalized = var(local_coverage_adaptive_simple_non_normalized, na.rm = TRUE),
  Adaptive_Simple_eGARCH_NonNormalized = var(local_coverage_adaptive_Esimple_non_normalized, na.rm = TRUE),
  Adaptive_Simple_GARCH_Normalized = var(local_coverage_adaptive_simple_normalized, na.rm = TRUE),
  Adaptive_Simple_eGARCH_Normalized = var(local_coverage_adaptive_Esimple_normalized, na.rm = TRUE),
  
  # Adaptive Momentum
  Adaptive_Momentum_GARCH = var(local_coverage_adaptive_momentum, na.rm = TRUE),
  Adaptive_Momentum_eGARCH = var(local_coverage_adaptive_Emomentum, na.rm = TRUE)
)

# Step 1: Best Model for Standard CP
standard_models <- c("Standard_CP_GARCH_NonNormalized", "Standard_CP_eGARCH_NonNormalized", 
                     "Standard_CP_GARCH_Normalized", "Standard_CP_eGARCH_Normalized")
standard_best <- standard_models[which.min(abs(unlist(model_means[standard_models]) - 0.9) + 
                                             unlist(model_variances[standard_models]))]

# Step 2: Best Model for Adaptive Simple
adaptive_simple_models <- c("Adaptive_Simple_GARCH_NonNormalized", "Adaptive_Simple_eGARCH_NonNormalized", 
                            "Adaptive_Simple_GARCH_Normalized", "Adaptive_Simple_eGARCH_Normalized")
adaptive_simple_best <- adaptive_simple_models[which.min(abs(unlist(model_means[adaptive_simple_models]) - 0.9) + 
                                                           unlist(model_variances[adaptive_simple_models]))]

# Step 3: Best Model for Adaptive Momentum
adaptive_momentum_models <- c("Adaptive_Momentum_GARCH", "Adaptive_Momentum_eGARCH")
adaptive_momentum_best <- adaptive_momentum_models[which.min(abs(unlist(model_means[adaptive_momentum_models]) - 0.9) + 
                                                               unlist(model_variances[adaptive_momentum_models]))]

# Step 4: Compare the Best Models from Each Category
final_models <- c(standard_best, adaptive_simple_best, adaptive_momentum_best)
final_best <- final_models[which.min(abs(unlist(model_means[final_models]) - 0.9) + 
                                       unlist(model_variances[final_models]))]

# Print Results
cat("Best model for Standard CP:", standard_best, "with mean:", model_means[[standard_best]], ", variance:", model_variances[[standard_best]], "\n")
## Best model for Standard CP: Standard_CP_GARCH_Normalized with mean: 0.9214204 , variance: 0.0001706544
cat("Best model for Adaptive Simple:", adaptive_simple_best, "with mean:", model_means[[adaptive_simple_best]], ", variance:", model_variances[[adaptive_simple_best]], "\n")
## Best model for Adaptive Simple: Adaptive_Simple_GARCH_Normalized with mean: 0.9025222 , variance: 3.857663e-05
cat("Best model for Adaptive Momentum:", adaptive_momentum_best, "with mean:", model_means[[adaptive_momentum_best]], ", variance:", model_variances[[adaptive_momentum_best]], "\n")
## Best model for Adaptive Momentum: Adaptive_Momentum_eGARCH with mean: 0.9024282 , variance: 4.937458e-05
cat("Best overall model for Gamma = 0.05:", final_best, "with mean:", model_means[[final_best]], ", variance:", model_variances[[final_best]], "\n")
## Best overall model for Gamma = 0.05: Adaptive_Momentum_eGARCH with mean: 0.9024282 , variance: 4.937458e-05

Results for Gamma = 0.05

For \(\gamma = 0.05\), the best-performing models for each strategy were identified based on the smallest absolute difference between the mean local coverage rate and the target (\(0.9\)), combined with the lowest variance.

Among the selected models, the Adaptive Momentum EGARCH strategy emerged as the best overall performer.


Next Steps

We will now evaluate the performance of the strategies for other values of \(\gamma\), specifically \(\gamma = 0.03\) and \(\gamma = 0.07\), to assess how the adaptation rate impacts the results and identify the optimal \(\gamma\) value.

# Analyze the impact of different gamma values on model selection

### Volatility prediction GAMMA= 0.03
# GARCH
GarchConformalForcasting_003 <- function(returns, alpha = 0.1, gamma = 0.03, lookback = 1250, 
                                         garchP = 1, garchQ = 1, startUp = 100, verbose = FALSE, 
                                         updateMethod, momentumBW = 0.8, scoreType ) {
  myT <- length(returns)
  T0 <- max(startUp, lookback)
  garchSpec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), 
                          variance.model = list(model = "sGARCH", garchOrder = c(garchP, garchQ)), 
                          distribution.model = "norm")
  alphat <- alpha
  ### Initialize data storage variables
  errSeqOC <- rep(0, myT - T0 + 1)
  errSeqNC <- rep(0, myT - T0 + 1)
  alphaSequence <- rep(alpha, myT - T0 + 1)
  scores <- rep(0, myT - T0 + 1)
  
  for (t in T0:myT) {
    if (verbose) {
      print(t)
    }
    ### Fit GARCH model and compute new conformity score
    garchFit <- ugarchfit(garchSpec, returns[(t - lookback + 1):(t - 1)], solver = "hybrid")
    sigmaNext <- sigma(ugarchforecast(garchFit, n.ahead = 1))
    if (scoreType == "Simple") {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2)
    } else {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2) / sigmaNext^2
    }
    
    recentScores <- scores[max(t - T0 + 1 - lookback + 1, 1):(t - T0)]
    
    ### Compute errt for both methods
    errSeqOC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alphat))))
    errSeqNC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alpha))))
    
    ### Update alphat
    alphaSequence[t - T0 + 1] <- alphat
    if (updateMethod == "Simple") {
      alphat <- alphat + gamma * (alpha - errSeqOC[t - T0 + 1])
    } else if (updateMethod == "Momentum") {
      w <- rev(momentumBW^(1:(t - T0 + 1)))
      w <- w / sum(w)
      alphat <- alphat + gamma * (alpha - sum(errSeqOC[1:(t - T0 + 1)] * w))
    }
    if (t %% 100 == 0) {
      print(sprintf("Done %g steps", t))
    }
  }
  
  return(list(alphaSequence = alphaSequence, errSeqOC = errSeqOC, errSeqNC = errSeqNC))
}

results_Studentized_003 <- GarchConformalForcasting_003(returns = returns, alpha = 0.1, gamma = 0.03,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
results_Simple_003 <- GarchConformalForcasting_003(returns = returns, alpha = 0.1, gamma = 0.03,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Simple")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
date_003 <- index(aapl_prices)[1250:length(index(aapl_prices))]
alphaSequence_St_003 <- results_Studentized_003[[1]]
errSeqOC_St_003 <- results_Studentized_003[[2]]; mean(errSeqOC_St_003)
## [1] 0.09753695
errSeqNC_St_003 <- results_Studentized_003[[3]]; mean(errSeqNC_St_003)
## [1] 0.08078818
alphaSequence_Si_003 <- results_Simple_003[[1]]
errSeqOC_Si_003 <- results_Simple_003[[2]]; mean(errSeqOC_Si_003)
## [1] 0.08866995
errSeqNC_Si_003 <- results_Simple_003[[3]]; mean(errSeqNC_Si_003)
## [1] 0.06699507
results_StudentizedM_003 <- GarchConformalForcasting_003(returns = returns, alpha = 0.1, gamma = 0.03,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Momentum",scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
alphaSequence_StM_003 <- results_StudentizedM_003[[1]]
errSeqOC_StM_003 <- results_StudentizedM_003[[2]]; mean(errSeqOC_StM_003)
## [1] 0.0955665
errSeqNC_StM_003 <- results_StudentizedM_003[[3]]; mean(errSeqNC_StM_003)
## [1] 0.08078818
plot(date_003, cummean(errSeqNC_St_003), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: GARCH Strategies Studentized (gamma = 0.03)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_003, cummean(errSeqOC_St_003), col = "blue")
lines(date_003, cummean(errSeqOC_StM_003), col = "red",type = "s")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple", "Adaptive Momentum"), 
       col = c("black", "blue", "red"), lty = 1)

plot(date_003, cummean(errSeqNC_Si_003), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: GARCH Strategies Simple (gamma = 0.03)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_003, cummean(errSeqOC_Si_003), col = "blue")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple"), 
       col = c("black", "blue"), lty = 1)

compute_local_coverage <- function(errors, window = 250) {
  1 - rollmean(errors, k = window, fill = NA)
}

window_size <- 250

# Generate Local Coverage Rates for Normalized and Non-Normalized Scores

# Local Coverage for Scores Non-Normalized
local_coverage_fixed_non_normalized_003 <- compute_local_coverage(errSeqNC_Si_003, window = window_size)
local_coverage_adaptive_simple_non_normalized_003 <- compute_local_coverage(errSeqOC_Si_003, window = window_size)

# Local Coverage for Scores Normalized
local_coverage_fixed_normalized_003 <- compute_local_coverage(errSeqNC_St_003, window = window_size)
local_coverage_adaptive_simple_normalized_003 <- compute_local_coverage(errSeqOC_St_003, window = window_size)
local_coverage_adaptive_momentum_003 <- compute_local_coverage(errSeqOC_StM_003, window = window_size)

# Plot Local Coverage Rates for Non-Normalized Scores
plot(date_003, local_coverage_adaptive_simple_non_normalized_003, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "GARCH Local Coverage Rates: Non-Normalized (gamma = 0.03)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_003, local_coverage_fixed_non_normalized_003, col = "red")
abline(h = 0.9, col = "black", lty = 2)
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha"), 
       col = c("blue", "red"), lty = 1)

# Plot Local Coverage Rates for Normalized Scores
plot(date_003, local_coverage_adaptive_simple_normalized_003, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "GARCH Local Coverage Rates: Normalized (gamma = 0.03)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_003, local_coverage_fixed_normalized_003, col = "red")
lines(date_003, local_coverage_adaptive_momentum_003, col = "green", type="s")
abline(h = 0.9, col = "black", lty = 2)  # Target coverage level
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha", "Adaptive Momentum"), 
       col = c("blue", "red", "green"), lty = 1)

Results for GARCH with Gamma = 0.03

Error Rates: Studentized

For studentized conformity scores, the Adaptive Simple and Adaptive Momentum strategies using the GARCH model closely approach the target error rate (\(\alpha = 0.1\)), with nearly identical performance. The Standard CP strategy, as observed previously, remains farther from the target throughout the period.


Error Rates: Non-Normalized

For non-normalized conformity scores, the Adaptive Simple strategy demonstrates significant improvement in convergence to the target error rate compared to the Standard CP, which performs worse consistently.


Local Coverage Rates: Non-Normalized

The local coverage rates for non-normalized conformity scores show a similar pattern to the previous cases. However, compared to \(\gamma = 0.05\), both strategies deviate more noticeably from the target level (\(90\%\)). The Standard CP remains significantly above the target, while the Adaptive Simple performs better but shows less alignment to the target compared to \(\gamma = 0.05\).


Local Coverage Rates: Normalized

For normalized conformity scores, the Adaptive Simple and Adaptive Momentum strategies using GARCH exhibit nearly identical behavior, oscillating closely around the target level (\(90\%\)). However, both strategies appear to deviate slightly more from the target compared to the results for \(\gamma = 0.05\). The Standard CP strategy remains consistently above the target.

### Volatility prediction GAMMA=0.03
# eGARCH
EgarchConformalForcasting_003 <- function(returns, alpha = 0.1, gamma = 0.03, lookback = 1250, 
                                          garchP = 1, garchQ = 1, startUp = 100, verbose = FALSE, 
                                          updateMethod, momentumBW = 0.8, scoreType ) {
  myT <- length(returns)
  T0 <- max(startUp, lookback)
  
  # Updated to eGARCH specification
  garchSpec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), 
                          variance.model = list(model = "eGARCH", garchOrder = c(garchP, garchQ)), 
                          distribution.model = "norm")
  
  alphat <- alpha
  errSeqOC <- rep(0, myT - T0 + 1)
  errSeqNC <- rep(0, myT - T0 + 1)
  alphaSequence <- rep(alpha, myT - T0 + 1)
  scores <- rep(0, myT - T0 + 1)
  
  for (t in T0:myT) {
    if (verbose) {
      print(t)
    }
    # Fit eGARCH model and compute conformity score
    garchFit <- ugarchfit(garchSpec, returns[(t - lookback + 1):(t - 1)], solver = "hybrid")
    sigmaNext <- sigma(ugarchforecast(garchFit, n.ahead = 1))
    
    if (scoreType == "Simple") {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2)
    } else {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2) / sigmaNext^2
    }
    
    recentScores <- scores[max(t - T0 + 1 - lookback + 1, 1):(t - T0)]
    errSeqOC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alphat))))
    errSeqNC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alpha))))
    
    alphaSequence[t - T0 + 1] <- alphat
    if (updateMethod == "Simple") {
      alphat <- alphat + gamma * (alpha - errSeqOC[t - T0 + 1])
    } else if (updateMethod == "Momentum") {
      w <- rev(momentumBW^(1:(t - T0 + 1)))
      w <- w / sum(w)
      alphat <- alphat + gamma * (alpha - sum(errSeqOC[1:(t - T0 + 1)] * w))
    }
    if (t %% 100 == 0) {
      print(sprintf("Done %g steps", t))
    }
  }
  
  return(list(alphaSequence = alphaSequence, errSeqOC = errSeqOC, errSeqNC = errSeqNC))
}

results_EStudentized_003 <- EgarchConformalForcasting_003(returns = returns, alpha = 0.1, gamma = 0.03,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
results_ESimple_003 <- EgarchConformalForcasting_003(returns = returns, alpha = 0.1, gamma = 0.03,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Simple")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
date_003 <- index(aapl_prices)[1250:length(index(aapl_prices))]
alphaSequence_ESt_003 <- results_EStudentized_003[[1]]
errSeqOC_ESt_003 <- results_EStudentized_003[[2]]; mean(errSeqOC_ESt_003)
## [1] 0.09753695
errSeqNC_ESt_003 <- results_EStudentized_003[[3]]; mean(errSeqNC_ESt_003)
## [1] 0.06699507
alphaSequence_ESi_003 <- results_ESimple_003[[1]]
errSeqOC_ESi_003 <- results_ESimple_003[[2]]; mean(errSeqOC_ESi_003)
## [1] 0.08866995
errSeqNC_ESi_003 <- results_ESimple_003[[3]]; mean(errSeqNC_ESi_003)
## [1] 0.06403941
results_EStudentizedM_003 <- EgarchConformalForcasting_003(returns = returns, alpha = 0.1, gamma = 0.03,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Momentum",scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
alphaSequence_EStM_003 <- results_EStudentizedM_003[[1]]
errSeqOC_EStM_003 <- results_EStudentizedM_003[[2]]; mean(errSeqOC_EStM_003)
## [1] 0.0955665
errSeqNC_EStM_003 <- results_EStudentizedM_003[[3]]; mean(errSeqNC_EStM_003)
## [1] 0.06699507
plot(date_003, cummean(errSeqNC_ESt_003), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: eGARCH Strategies Studentized (gamma = 0.03)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_003, cummean(errSeqOC_ESt_003), col = "blue")
lines(date_003, cummean(errSeqOC_EStM_003), col = "red",type = "s")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple", "Adaptive Momentum"), 
       col = c("black", "blue", "red"), lty = 1)

plot(date_003, cummean(errSeqNC_ESi_003), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: eGARCH Strategies Simple (gamma = 0.03)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_003, cummean(errSeqOC_ESi_003), col = "blue")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple"), 
       col = c("black", "blue"), lty = 1)

compute_local_coverage <- function(errors, window = 250) {
  1 - rollmean(errors, k = window, fill = NA)
}

window_size <- 250
# Generate Local Coverage Rates for Normalized and Non-Normalized Scores

# Local Coverage for Scores Non-Normalized (eGARCH)
local_coverage_Efixed_non_normalized_003 <- compute_local_coverage(errSeqNC_ESi_003, window = window_size)
local_coverage_adaptive_Esimple_non_normalized_003 <- compute_local_coverage(errSeqOC_ESi_003, window = window_size)

# Local Coverage for Scores Normalized (eGARCH)
local_coverage_Efixed_normalized_003 <- compute_local_coverage(errSeqNC_ESt_003, window = window_size)
local_coverage_adaptive_Esimple_normalized_003 <- compute_local_coverage(errSeqOC_ESt_003, window = window_size)
local_coverage_adaptive_Emomentum_003 <- compute_local_coverage(errSeqOC_EStM_003, window = window_size)

# Plot Local Coverage Rates for Non-Normalized Scores (eGARCH)
plot(date_003, local_coverage_adaptive_Esimple_non_normalized_003, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "eGARCH Local Coverage Rates: Non-Normalized (gamma = 0.03)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_003, local_coverage_Efixed_non_normalized_003, col = "red")
abline(h = 0.9, col = "black", lty = 2)
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha"), 
       col = c("blue", "red"), lty = 1)

# Plot Local Coverage Rates for Normalized Scores (eGARCH)
plot(date_003, local_coverage_adaptive_Esimple_normalized_003, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "eGARCH Local Coverage Rates: Normalized (gamma = 0.03)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_003, local_coverage_Efixed_normalized_003, col = "red")
lines(date_003, local_coverage_adaptive_Emomentum_003, col = "green", type="s")
abline(h = 0.9, col = "black", lty = 2)  # Target coverage level
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha", "Adaptive Momentum"), 
       col = c("blue", "red", "green"), lty = 1)

Results for eGARCH with Gamma = 0.03

Error Rates: Studentized

For studentized conformity scores, the results obtained with the eGARCH model are very similar to those observed for GARCH. Both the Adaptive Simple and Adaptive Momentum strategies closely approach the target error rate (\(\alpha = 0.1\)), while the Standard CP strategy remains farther from the target throughout the period.


Error Rates: Non-Normalized

For non-normalized conformity scores, the Adaptive Simple strategy with eGARCH also demonstrates better convergence to the target error rate compared to the Standard CP, which performs worse consistently.


Local Coverage Rates: Non-Normalized

The local coverage rates for non-normalized conformity scores using eGARCH follow a similar pattern to those observed with GARCH. Compared to \(\gamma = 0.05\), both the Adaptive Simple and Standard CP strategies deviate more noticeably from the target level (\(90\%\)). The Standard CP remains significantly above the target, while the Adaptive Simple performs better but is slightly less aligned to the target compared to \(\gamma = 0.05\).


Local Coverage Rates: Normalized

For normalized conformity scores, the Adaptive Simple and Adaptive Momentum strategies using eGARCH show nearly identical performance, oscillating closely around the target level (\(90\%\)). However, both strategies deviate slightly more from the target compared to \(\gamma = 0.05\). The Standard CP strategy remains consistently above the target.

# Compare and Select the Best Models for Gamma = 0.03
model_means_003 <- list(
  # Standard CP
  Standard_CP_GARCH_Normalized = mean(local_coverage_fixed_normalized_003, na.rm = TRUE),
  
  # Adaptive Simple
  Adaptive_Simple_GARCH_Normalized = mean(local_coverage_adaptive_simple_normalized_003, na.rm = TRUE),
  
  # Adaptive Momentum
  Adaptive_Momentum_eGARCH = mean(local_coverage_adaptive_Emomentum_003, na.rm = TRUE)
)

# Compare models using both mean and variance of local coverage rates
model_variances_003 <- list(
  # Standard CP
  Standard_CP_GARCH_Normalized = var(local_coverage_fixed_normalized_003, na.rm = TRUE),
  
  # Adaptive Simple
  Adaptive_Simple_GARCH_Normalized = var(local_coverage_adaptive_simple_normalized_003, na.rm = TRUE),
  
  # Adaptive Momentum
  Adaptive_Momentum_eGARCH = var(local_coverage_adaptive_Emomentum_003, na.rm = TRUE)
)

# Step 1: Best Model for Gamma = 0.03
models_003 <- c("Standard_CP_GARCH_Normalized", 
                "Adaptive_Simple_GARCH_Normalized", 
                "Adaptive_Momentum_eGARCH")

best_model_003 <- models_003[which.min(
  abs(unlist(model_means_003[models_003]) - 0.9) + 
    unlist(model_variances_003[models_003])
)]

# Print Results for Gamma = 0.03
cat("Best model for Gamma = 0.03:", best_model_003, "with mean:", 
    model_means_003[[best_model_003]], ", variance:", 
    model_variances_003[[best_model_003]], "\n")
## Best model for Gamma = 0.03: Adaptive_Simple_GARCH_Normalized with mean: 0.9031802 , variance: 8.091652e-05

Results for Gamma = 0.03: Model Comparison

After comparing the results of the three best models identified for \(\gamma = 0.05\), the Adaptive_Simple_GARCH_Normalized model was selected as the best-performing strategy for \(\gamma = 0.03\).

### Volatility prediction GAMMA= 0.07
# GARCH
GarchConformalForcasting_007 <- function(returns, alpha = 0.1, gamma = 0.07, lookback = 1250, 
                                         garchP = 1, garchQ = 1, startUp = 100, verbose = FALSE, 
                                         updateMethod, momentumBW = 0.8, scoreType ) {
  myT <- length(returns)
  T0 <- max(startUp, lookback)
  garchSpec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), 
                          variance.model = list(model = "sGARCH", garchOrder = c(garchP, garchQ)), 
                          distribution.model = "norm")
  alphat <- alpha
  ### Initialize data storage variables
  errSeqOC <- rep(0, myT - T0 + 1)
  errSeqNC <- rep(0, myT - T0 + 1)
  alphaSequence <- rep(alpha, myT - T0 + 1)
  scores <- rep(0, myT - T0 + 1)
  
  for (t in T0:myT) {
    if (verbose) {
      print(t)
    }
    ### Fit GARCH model and compute new conformity score
    garchFit <- ugarchfit(garchSpec, returns[(t - lookback + 1):(t - 1)], solver = "hybrid")
    sigmaNext <- sigma(ugarchforecast(garchFit, n.ahead = 1))
    if (scoreType == "Simple") {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2)
    } else {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2) / sigmaNext^2
    }
    
    recentScores <- scores[max(t - T0 + 1 - lookback + 1, 1):(t - T0)]
    
    ### Compute errt for both methods
    errSeqOC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alphat))))
    errSeqNC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alpha))))
    
    ### Update alphat
    alphaSequence[t - T0 + 1] <- alphat
    if (updateMethod == "Simple") {
      alphat <- alphat + gamma * (alpha - errSeqOC[t - T0 + 1])
    } else if (updateMethod == "Momentum") {
      w <- rev(momentumBW^(1:(t - T0 + 1)))
      w <- w / sum(w)
      alphat <- alphat + gamma * (alpha - sum(errSeqOC[1:(t - T0 + 1)] * w))
    }
    if (t %% 100 == 0) {
      print(sprintf("Done %g steps", t))
    }
  }
  
  return(list(alphaSequence = alphaSequence, errSeqOC = errSeqOC, errSeqNC = errSeqNC))
}

results_Studentized_007 <- GarchConformalForcasting_007(returns = returns, alpha = 0.1, gamma = 0.07,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
results_Simple_007 <- GarchConformalForcasting_007(returns = returns, alpha = 0.1, gamma = 0.07,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Simple")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
date_007 <- index(aapl_prices)[1250:length(index(aapl_prices))]
alphaSequence_St_007 <- results_Studentized_007[[1]]
errSeqOC_St_007 <- results_Studentized_007[[2]]; mean(errSeqOC_St_007)
## [1] 0.09852217
errSeqNC_St_007 <- results_Studentized_007[[3]]; mean(errSeqNC_St_007)
## [1] 0.08078818
alphaSequence_Si_007 <- results_Simple_007[[1]]
errSeqOC_Si_007 <- results_Simple_007[[2]]; mean(errSeqOC_Si_007)
## [1] 0.09261084
errSeqNC_Si_007 <- results_Simple_007[[3]]; mean(errSeqNC_Si_007)
## [1] 0.06699507
results_StudentizedM_007 <- GarchConformalForcasting_007(returns = returns, alpha = 0.1, gamma = 0.07,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Momentum",scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
alphaSequence_StM_007 <- results_StudentizedM_007[[1]]
errSeqOC_StM_007 <- results_StudentizedM_007[[2]]; mean(errSeqOC_StM_007)
## [1] 0.09852217
errSeqNC_StM_007 <- results_StudentizedM_007[[3]]; mean(errSeqNC_StM_007)
## [1] 0.08078818
plot(date_007, cummean(errSeqNC_St_007), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: GARCH Strategies Studentized (gamma = 0.07)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_007, cummean(errSeqOC_St_007), col = "blue")
lines(date_007, cummean(errSeqOC_StM_007), col = "red",type = "s")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple", "Adaptive Momentum"), 
       col = c("black", "blue", "red"), lty = 1)

plot(date_007, cummean(errSeqNC_Si_007), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: GARCH Strategies Simple (gamma = 0.07)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_007, cummean(errSeqOC_Si_007), col = "blue")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple"), 
       col = c("black", "blue"), lty = 1)

compute_local_coverage <- function(errors, window = 250) {
  1 - rollmean(errors, k = window, fill = NA)
}

window_size <- 250

# Generate Local Coverage Rates for Normalized and Non-Normalized Scores

# Local Coverage for Scores Non-Normalized
local_coverage_fixed_non_normalized_007 <- compute_local_coverage(errSeqNC_Si_007, window = window_size)
local_coverage_adaptive_simple_non_normalized_007 <- compute_local_coverage(errSeqOC_Si_007, window = window_size)

# Local Coverage for Scores Normalized
local_coverage_fixed_normalized_007 <- compute_local_coverage(errSeqNC_St_007, window = window_size)
local_coverage_adaptive_simple_normalized_007 <- compute_local_coverage(errSeqOC_St_007, window = window_size)
local_coverage_adaptive_momentum_007 <- compute_local_coverage(errSeqOC_StM_007, window = window_size)

# Plot Local Coverage Rates for Non-Normalized Scores
plot(date_007, local_coverage_adaptive_simple_non_normalized_007, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "GARCH Local Coverage Rates: Non-Normalized (gamma = 0.07)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_007, local_coverage_fixed_non_normalized_007, col = "red")
abline(h = 0.9, col = "black", lty = 2)
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha"), 
       col = c("blue", "red"), lty = 1)

# Plot Local Coverage Rates for Normalized Scores
plot(date_007, local_coverage_adaptive_simple_normalized_007, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "GARCH Local Coverage Rates: Normalized (gamma = 0.07)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_007, local_coverage_fixed_normalized_007, col = "red")
lines(date_007, local_coverage_adaptive_momentum_007, col = "green", type="s")
abline(h = 0.9, col = "black", lty = 2)  # Target coverage level
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha", "Adaptive Momentum"), 
       col = c("blue", "red", "green"), lty = 1)

Results for GARCH with Gamma = 0.07

Error Rates: Studentized

For studentized conformity scores, the results with \(\gamma = 0.07\) closely resemble those observed for \(\gamma = 0.05\). Both the Adaptive Simple and Adaptive Momentum strategies perform very well, approaching the target error rate (\(\alpha = 0.1\)) with nearly identical behavior. The Standard CP strategy, as usual, remains farther from the target.


Error Rates: Non-Normalized

For non-normalized conformity scores, the Adaptive Simple strategy demonstrates clear superiority in approaching the target error rate compared to the Standard CP, which consistently shows larger deviations from the target.


Local Coverage Rates: Non-Normalized

The local coverage rates for non-normalized conformity scores using \(\gamma = 0.07\) align closely with the results for \(\gamma = 0.05\). Both strategies oscillate around the target level (\(90\%\)), but the Adaptive Simple strategy is better aligned compared to the Standard CP, which consistently remains above the target.


Local Coverage Rates: Normalized

For normalized conformity scores, the Adaptive Simple and Adaptive Momentum strategies demonstrate strong performance, oscillating closely around the target coverage level (\(90\%\)), and their results are nearly indistinguishable. The Standard CP strategy, as in previous cases, remains consistently above the target.

### Volatility prediction GAMMA= 0.07
# eGARCH
EgarchConformalForcasting_007 <- function(returns, alpha = 0.1, gamma = 0.07, lookback = 1250, 
                                          garchP = 1, garchQ = 1, startUp = 100, verbose = FALSE, 
                                          updateMethod, momentumBW = 0.8, scoreType ) {
  myT <- length(returns)
  T0 <- max(startUp, lookback)
  
  # Updated to eGARCH specification
  garchSpec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean = FALSE), 
                          variance.model = list(model = "eGARCH", garchOrder = c(garchP, garchQ)), 
                          distribution.model = "norm")
  
  alphat <- alpha
  errSeqOC <- rep(0, myT - T0 + 1)
  errSeqNC <- rep(0, myT - T0 + 1)
  alphaSequence <- rep(alpha, myT - T0 + 1)
  scores <- rep(0, myT - T0 + 1)
  
  for (t in T0:myT) {
    if (verbose) {
      print(t)
    }
    # Fit eGARCH model and compute conformity score
    garchFit <- ugarchfit(garchSpec, returns[(t - lookback + 1):(t - 1)], solver = "hybrid")
    sigmaNext <- sigma(ugarchforecast(garchFit, n.ahead = 1))
    
    if (scoreType == "Simple") {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2)
    } else {
      scores[t - T0 + 1] <- abs(returns[t]^2 - sigmaNext^2) / sigmaNext^2
    }
    
    recentScores <- scores[max(t - T0 + 1 - lookback + 1, 1):(t - T0)]
    errSeqOC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alphat))))
    errSeqNC[t - T0 + 1] <- as.numeric(scores[t - T0 + 1] > quantile(recentScores, probs = max(0, min(1, 1 - alpha))))
    
    alphaSequence[t - T0 + 1] <- alphat
    if (updateMethod == "Simple") {
      alphat <- alphat + gamma * (alpha - errSeqOC[t - T0 + 1])
    } else if (updateMethod == "Momentum") {
      w <- rev(momentumBW^(1:(t - T0 + 1)))
      w <- w / sum(w)
      alphat <- alphat + gamma * (alpha - sum(errSeqOC[1:(t - T0 + 1)] * w))
    }
    if (t %% 100 == 0) {
      print(sprintf("Done %g steps", t))
    }
  }
  
  return(list(alphaSequence = alphaSequence, errSeqOC = errSeqOC, errSeqNC = errSeqNC))
}

results_EStudentized_007 <- EgarchConformalForcasting_007(returns = returns, alpha = 0.1, gamma = 0.07,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
results_ESimple_007 <- EgarchConformalForcasting_007(returns = returns, alpha = 0.1, gamma = 0.07,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Simple", scoreType ="Simple")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
date_007 <- index(aapl_prices)[1250:length(index(aapl_prices))]
alphaSequence_ESt_007 <- results_EStudentized_007[[1]]
errSeqOC_ESt_007 <- results_EStudentized_007[[2]]; mean(errSeqOC_ESt_007)
## [1] 0.09852217
errSeqNC_ESt_007 <- results_EStudentized_007[[3]]; mean(errSeqNC_ESt_007)
## [1] 0.06699507
alphaSequence_ESi_007 <- results_ESimple_007[[1]]
errSeqOC_ESi_007 <- results_ESimple_007[[2]]; mean(errSeqOC_ESi_007)
## [1] 0.09359606
errSeqNC_ESi_007 <- results_ESimple_007[[3]]; mean(errSeqNC_ESi_007)
## [1] 0.06403941
results_EStudentizedM_007 <- EgarchConformalForcasting_007(returns = returns, alpha = 0.1, gamma = 0.07,lookback=1250,garchP=1,garchQ=1,startUp = 100,verbose=FALSE,updateMethod="Momentum",scoreType ="Studentized")
## [1] "Done 1300 steps"
## [1] "Done 1400 steps"
## [1] "Done 1500 steps"
## [1] "Done 1600 steps"
## [1] "Done 1700 steps"
## [1] "Done 1800 steps"
## [1] "Done 1900 steps"
## [1] "Done 2000 steps"
## [1] "Done 2100 steps"
## [1] "Done 2200 steps"
alphaSequence_EStM_007 <- results_EStudentizedM_007[[1]]
errSeqOC_EStM_007 <- results_EStudentizedM_007[[2]]; mean(errSeqOC_EStM_007)
## [1] 0.09852217
errSeqNC_EStM_007 <- results_EStudentizedM_007[[3]]; mean(errSeqNC_EStM_007)
## [1] 0.06699507
plot(date_007, cummean(errSeqNC_ESt_007), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: eGARCH Strategies Studentized (gamma = 0.07)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_007, cummean(errSeqOC_ESt_007), col = "blue")
lines(date_007, cummean(errSeqOC_EStM_007), col = "red",type = "s")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple", "Adaptive Momentum"), 
       col = c("black", "blue", "red"), lty = 1)

plot(date_007, cummean(errSeqNC_ESi_007), type = "l", col = "black", ylim = c(0, 0.5), 
     main = "Error Rates: eGARCH Strategies Simple (gamma = 0.07)", xlab = "Date", ylab = "Cumulative Error Rate")
lines(date_007, cummean(errSeqOC_ESi_007), col = "blue")
abline(h = 0.1, col = "green", lty = 2)
legend("topright", legend = c("Standard CP", "Adaptive Simple"), 
       col = c("black", "blue"), lty = 1)

compute_local_coverage <- function(errors, window = 250) {
  1 - rollmean(errors, k = window, fill = NA)
}

window_size <- 250
# Generate Local Coverage Rates for Normalized and Non-Normalized Scores

# Local Coverage for Scores Non-Normalized (eGARCH)
local_coverage_Efixed_non_normalized_007 <- compute_local_coverage(errSeqNC_ESi_007, window = window_size)
local_coverage_adaptive_Esimple_non_normalized_007 <- compute_local_coverage(errSeqOC_ESi_007, window = window_size)

# Local Coverage for Scores Normalized (eGARCH)
local_coverage_Efixed_normalized_007 <- compute_local_coverage(errSeqNC_ESt_007, window = window_size)
local_coverage_adaptive_Esimple_normalized_007 <- compute_local_coverage(errSeqOC_ESt_007, window = window_size)
local_coverage_adaptive_Emomentum_007 <- compute_local_coverage(errSeqOC_EStM_007, window = window_size)

# Plot Local Coverage Rates for Non-Normalized Scores (eGARCH)
plot(date_007, local_coverage_adaptive_Esimple_non_normalized_007, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "eGARCH Local Coverage Rates: Non-Normalized (gamma = 0.07)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_007, local_coverage_Efixed_non_normalized_007, col = "red")
abline(h = 0.9, col = "black", lty = 2)
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha"), 
       col = c("blue", "red"), lty = 1)

# Plot Local Coverage Rates for Normalized Scores (eGARCH)
plot(date_007, local_coverage_adaptive_Esimple_normalized_007, type = "l", col = "blue", ylim = c(0.8, 1),
     main = "eGARCH Local Coverage Rates: Normalized (gamma = 0.07)", xlab = "Date", ylab = "Local Coverage Level")
lines(date_007, local_coverage_Efixed_normalized_007, col = "red")
lines(date_007, local_coverage_adaptive_Emomentum_007, col = "green", type="s")
abline(h = 0.9, col = "black", lty = 2)  # Target coverage level
legend("bottomright", legend = c("Adaptive Simple", "Fixed Alpha", "Adaptive Momentum"), 
       col = c("blue", "red", "green"), lty = 1)

Results for eGARCH with Gamma = 0.07

Error Rates: Studentized

For studentized conformity scores, the results with \(\gamma = 0.07\) for eGARCH closely align with those observed for GARCH. Both the Adaptive Simple and Adaptive Momentum strategies approach the target error rate (\(\alpha = 0.1\)) effectively, with nearly identical performance. The Standard CP strategy, as before, remains farther from the target.


Error Rates: Non-Normalized

For non-normalized conformity scores, the Adaptive Simple strategy with eGARCH shows clear superiority in converging to the target error rate compared to the Standard CP, which consistently deviates more significantly.


Local Coverage Rates: Non-Normalized

The local coverage rates for non-normalized conformity scores using eGARCH are highly similar to those seen with GARCH. Both the Adaptive Simple and Standard CP strategies oscillate around the target level (\(90\%\)), but the Adaptive Simple strategy is better aligned, while the Standard CP remains consistently above the target.


Local Coverage Rates: Normalized

For normalized conformity scores, the Adaptive Simple and Adaptive Momentum strategies using eGARCH perform nearly identically, oscillating closely around the target coverage level (\(90\%\)). The Standard CP strategy, as in previous cases, remains consistently above the target.

# Compare and Select the Best Models for Gamma = 0.07
model_means_007 <- list(
  # Standard CP
  Standard_CP_GARCH_Normalized = mean(local_coverage_fixed_normalized_007, na.rm = TRUE),
  
  # Adaptive Simple
  Adaptive_Simple_GARCH_Normalized = mean(local_coverage_adaptive_simple_normalized_007, na.rm = TRUE),
  
  # Adaptive Momentum
  Adaptive_Momentum_eGARCH = mean(local_coverage_adaptive_Emomentum_007, na.rm = TRUE)
)

# Compare models using both mean and variance of local coverage rates
model_variances_007 <- list(
  # Standard CP
  Standard_CP_GARCH_Normalized = var(local_coverage_fixed_normalized_007, na.rm = TRUE),
  
  # Adaptive Simple
  Adaptive_Simple_GARCH_Normalized = var(local_coverage_adaptive_simple_normalized_007, na.rm = TRUE),
  
  # Adaptive Momentum
  Adaptive_Momentum_eGARCH = var(local_coverage_adaptive_Emomentum_007, na.rm = TRUE)
)

# Step 1: Best Model for Gamma = 0.07
models_007 <- c("Standard_CP_GARCH_Normalized", 
                "Adaptive_Simple_GARCH_Normalized", 
                "Adaptive_Momentum_eGARCH")

best_model_007 <- models_007[which.min(
  abs(unlist(model_means_007[models_007]) - 0.9) + 
    unlist(model_variances_007[models_007])
)]

# Print Results for Gamma = 0.07
cat("Best model for Gamma = 0.07:", best_model_007, "with mean:", 
    model_means_007[[best_model_007]], ", variance:", 
    model_variances_007[[best_model_007]], "\n")
## Best model for Gamma = 0.07: Adaptive_Momentum_eGARCH with mean: 0.9015614 , variance: 4.062303e-05

Results for Gamma = 0.07: Model Comparison

After comparing the results of the three best models identified for \(\gamma = 0.05\), the Adaptive_Momentum_eGARCH model was selected as the best-performing strategy for \(\gamma = 0.07\).

# Compare Best Models Across All Gamma Values
final_model_means <- list(
  Gamma_0.03 = 0.9031802, # Adaptive_Simple_GARCH_Normalized for Gamma = 0.03
  Gamma_0.05 = 0.9024282, # Adaptive_Momentum_eGARCH for Gamma = 0.05
  Gamma_0.07 = 0.9015614  # Adaptive_Momentum_eGARCH for Gamma = 0.07
)

final_model_variances <- list(
  Gamma_0.03 = 8.091652e-05, # Variance for Gamma = 0.03
  Gamma_0.05 = 4.937458e-05, # Variance for Gamma = 0.05
  Gamma_0.07 = 4.062303e-05  # Variance for Gamma = 0.07
)

# Determine the Best Overall Model
final_models <- c("Gamma_0.03", "Gamma_0.05", "Gamma_0.07")

best_overall_model <- final_models[which.min(
  abs(unlist(final_model_means) - 0.9) + 
    unlist(final_model_variances)
)]

# Print the Best Overall Model
cat("Best overall model across all gamma values:", best_overall_model, "with mean:", 
    final_model_means[[best_overall_model]], ", variance:", 
    final_model_variances[[best_overall_model]], "\n")
## Best overall model across all gamma values: Gamma_0.07 with mean: 0.9015614 , variance: 4.062303e-05

Final Conclusion: Best Model Across All Gamma Values

The Adaptive_Momentum_eGARCH model with \(\gamma = 0.07\) was the best-performing model overall based on the criteria we used.

The method we used to select the best model (minimizing the difference between the mean local coverage rate and the target, along with the variance of the local coverage rate) was chosen because we considered it representative of the best-performing model. However, we recognize that this approach may be incomplete and could be improved.

Part 2 - Point a

Formal Description of AgACI

Background on Adaptive Conformal Inference (ACI)

The Adaptive Conformal Inference (ACI) framework is designed to create valid prediction intervals for time series data. The key challenge lies in dynamically adapting these intervals to maintain the desired coverage level \(1 - \alpha\), where \(\alpha\) is the miscoverage rate.

In ACI, this adaptation is controlled by an adaptation rate parameter, \(\gamma > 0\), which determines how quickly the intervals respond to changes in the underlying data distribution. Selecting the correct value of \(\gamma\) is crucial, as:

  1. Over-adaptation: Intervals become too sensitive to noise, leading to instability.
  2. Under-adaptation: Intervals fail to adjust to distributional shifts, causing poor coverage.

Motivation for Aggregated Adaptive Conformal Inference (AgACI)

The Aggregated Adaptive Conformal Inference (AgACI) method addresses the issue of selecting a single \(\gamma\) by introducing a parameter-free ensemble approach. Instead of relying on a fixed \(\gamma\), AgACI combines predictions from multiple experts, each using a distinct \(\gamma_k\).

This adaptive aggregation ensures robust performance across various conditions, mitigating the limitations of ACI.

Components of AgACI

Experts with Different \(\gamma\) Values

  • AgACI defines a set of \(K\) experts, each using a unique \(\gamma_k, k = 1, \dots, K\).
  • These \(\gamma_k\) values span a range of possible adaptation rates, from highly responsive to highly stable.

Expert Weights

  • At each time step \(t\), each expert is assigned a weight \(w_k(t)\) representing its confidence level.
  • Initially, weights are uniformly distributed across all experts: \(w_k(0) = \frac{1}{K}\).

Prediction Interval Aggregation

  • Each expert \(k\) produces a prediction interval \([L_k(t), U_k(t)]\) based on its \(\gamma_k\).
  • The aggregated interval is computed as:

\[ [L(t), U(t)] = \left[ \sum_{k=1}^K w_k(t) L_k(t), \sum_{k=1}^K w_k(t) U_k(t) \right]. \]

Weight Updates

  • Weights are updated dynamically at each time step based on the experts’ past performance.
  • The pinball loss function is used to evaluate expert performance. For a quantile level \(\tau\):

\[ \mathcal{L}(y_t, \hat{y}_t) = \begin{cases} \tau (y_t - \hat{y}_t), & y_t > \hat{y}_t, \\ (1 - \tau)(\hat{y}_t - y_t), & y_t \leq \hat{y}_t. \end{cases} \]

  • Experts with better performance (lower loss) are assigned higher weights in subsequent time steps.

#Point b

Using EGARCH in AgACI

EGARCH (Exponential Generalized Autoregressive Conditional Heteroskedasticity) models are specifically designed to capture volatility dynamics in financial time series.

Forecasting Volatility

EGARCH models the conditional variance of returns, capturing key characteristics such as volatility clustering and asymmetry. The forecasted variance, \(\hat{\sigma}_t^2\), is the central estimate of volatility.

Constructing Prediction Intervals

Prediction intervals can be defined using the forecasted variance scaled by quantiles of the residual distribution:

\[ \text{Interval} = \left[ \hat{\sigma}_t^2 \cdot Q_{\alpha/2}, \hat{\sigma}_t^2 \cdot Q_{1-\alpha/2} \right] \]

where \(Q_\alpha\) represents the quantile of the residuals.

Conformity Scores

Conformity scores measure how well the predictions align with observations. For EGARCH, the conformity score is given by:

\[ S_t = \frac{|r_t^2 - \hat{\sigma}_t^2|}{\hat{\sigma}_t^2} \]

Adaptive Updates

AgACI dynamically adjusts the size of prediction intervals based on past conformity scores, ensuring robust coverage in changing market conditions.

Advantages and Limitations of EGARCH

  • Advantages: EGARCH is designed for volatility modeling, capturing market dynamics such as asymmetry and clustering, and provides probabilistic forecasts.

  • Limitations: It assumes specific distributions (e.g., normal, Student’s t) and may struggle with extreme tail behavior compared to more flexible models.

Using Quantile Regression in AgACI

Quantile Regression offers a flexible and non-parametric alternative by modeling specific quantiles of the return distribution directly.

Modeling Return Quantiles

Quantile Regression predicts the lower and upper quantiles of future returns, forming the basis of prediction intervals:

\[ \text{Interval} = \left[ \hat{Q}_{\alpha/2}, \hat{Q}_{1-\alpha/2} \right] \]

Adaptive Updates

AgACI adjusts these quantiles based on observed conformity scores, enabling adaptive prediction intervals.

Advantages and Limitations of Quantile Regression

  • Advantages: Quantile Regression is flexible, robust to outliers, and well-suited for modeling extreme tail behavior.

  • Limitations: It does not explicitly account for volatility dynamics like clustering or asymmetry and focuses solely on return quantiles rather than variance.

Conclusion

The choice between EGARCH and Quantile Regression in AgACI depends on the specific objectives:

  • EGARCH excels in capturing volatility dynamics and market behaviors, making it ideal for volatility forecasting.

  • Quantile Regression provides flexibility and robustness for modeling extreme returns, making it a strong candidate for risk management.

Both approaches complement each other, offering valuable insights into financial time series depending on the context of the analysis.

Introduction

The Bernstein Online Aggregation (BOA) is an advanced algorithm for expert aggregation that dynamically combines the predictions from multiple models (or experts) to minimize overall prediction error. It is particularly well-suited for online learning scenarios where data arrives sequentially, and the performance of individual experts can vary over time. BOA is designed to be both adaptive and robust, ensuring that the aggregation process remains stable while quickly responding to changes in expert performance.

Theoretical Foundation

BOA is grounded in statistical learning theory and relies on Bernstein inequalities to control the trade-off between bias and variance in the aggregation process. Traditional aggregation methods may suffer from instability due to abrupt weight changes, but BOA addresses this issue by incorporating concentration inequalities that allow for smoother and more controlled weight updates.

Bernstein inequalities provide tighter control over the deviation between the true loss and the empirical loss, especially when dealing with bounded losses. This makes BOA particularly suitable for applications where prediction errors must be tightly managed.

Working Mechanism

Step-by-step process of BOA:

Initialization: - Each expert is assigned an initial weight, usually uniform across all experts (e.g., equal weights if no prior information is available).

Prediction: - At each time step \(t\), each expert provides a prediction for the next data point. - BOA computes a weighted average of these predictions using the current weights.

Loss Computation: - Once the true value is observed, BOA calculates the loss (error) for each expert. - The most common loss functions used are the squared loss for regression or the absolute loss for more robust applications.

Weight Update: - BOA updates the weights of each expert using a rule derived from Bernstein inequalities. - Experts that perform poorly receive exponentially decreasing weights, while well-performing experts are rewarded. - This update is smoother than in traditional algorithms, preventing extreme changes.

Aggregation: - The aggregated prediction is a weighted sum of expert predictions, where the weights reflect each expert’s historical performance.

Iteration: - This process repeats for each new data point, allowing BOA to adapt in real-time.

Mathematical Formulation

The weight update rule in BOA is designed to minimize the cumulative regret. The weight \(w_i(t)\) for expert \(i\) at time \(t\) is updated based on its cumulative loss using a function that penalizes high variance in predictions.

The general form of the weight update is:

\[ w_i(t+1) = \frac{w_i(t) \cdot \exp(-\eta \cdot L_i(t))}{Z(t)} \]

Where: - \(w_i(t)\) is the weight of expert \(i\) at time \(t\). - \(L_i(t)\) is the loss of expert \(i\) at time \(t\). - \(\eta\) is the learning rate, controlling how quickly the algorithm adapts. - \(Z(t)\) is a normalization factor ensuring that the weights sum to 1.

Advantages of BOA

  • Adaptive Learning: Quickly adapts to changes in expert performance, making it suitable for non-stationary environments like financial markets.
  • Robustness: Controls weight updates with Bernstein inequalities to avoid overreacting to short-term fluctuations.
  • Stability: Smooth weight updates prevent abrupt shifts that can lead to instability.
  • Regret Minimization: Designed to minimize regret, ensuring competitiveness with the best expert or combination of experts in hindsight.

Conclusion

The Bernstein Online Aggregation (BOA) algorithm is a powerful and adaptive method for combining expert predictions in online learning settings. Its integration of Bernstein inequalities allows for smooth and robust weight updates, making it especially effective in volatile and dynamic environments like financial forecasting. BOA provides a practical and efficient tool for expert aggregation, balancing performance and risk in prediction models.

library(opera)

agaciVolatilityForecasting <- function(returns, alpha = 0.1, 
                                       gammaGrid = c(0.03, 0.05, 0.07), 
                                       lookback = 1250, garchP = 1, garchQ = 1, 
                                       startUp = 100, verbose = FALSE, 
                                       updateMethod = "Momentum", momentumBW = 0.80) {
  myT <- length(returns)
  T0 <- max(startUp, lookback)
  
  
  garchSpec <- ugarchspec(
    mean.model = list(armaOrder = c(0, 0), include.mean = FALSE),
    variance.model = list(model = "eGARCH", garchOrder = c(garchP, garchQ)),
    distribution.model = "norm"  
  )
  
  # Initialize storage variables
  K <- length(gammaGrid)  # experts
  Tlen <- myT - T0
  
  
  garchForecastVec <- rep(NA, Tlen)
  
  scores          <- array(NA, dim = Tlen)
  alphaSeq        <- matrix(alpha, nrow = K, ncol = Tlen)
  
  alphaSeq2       <- rep(alpha, Tlen)  
  
  errSeqOC        <- matrix(NA, nrow = K, ncol = Tlen)  
  aggregatedalpha <- rep(NA, Tlen)
  err             <- rep(NA, Tlen)  
  
  lower_limit <- 0.001
  upper_limit <- 0.999
  
  
  
  for (t in T0:(myT - 1)) {
    
    # Fit GARCH 
    garchFit     <- ugarchfit(garchSpec, returns[(t - lookback + 1):(t - 1)], solver = "hybrid")
    sigmaNext    <- sigma(ugarchforecast(garchFit, n.ahead = 1))
    garchForecast <- sigmaNext^2  
    
    colIdx <- t - T0 + 1   
    
    
    garchForecastVec[colIdx] <- garchForecast
    
    # Conformity score
    scores[colIdx] <- abs(returns[t]^2 - sigmaNext^2) / sigmaNext^2  
    
    
    historical <- scores[max(1, colIdx - lookback + 1):(colIdx - 1)]
    
    
    for (g in seq_along(gammaGrid)) {
      gamma <- gammaGrid[g]
      
      if (t == T0) {
        thr <- quantile(historical, 1 - alpha)
        errSeqOC[g, colIdx] <- as.numeric(scores[colIdx] > thr)
        
        
      } else {
        thr <- quantile(historical, 1 - alphaSeq[g, colIdx - 1])
        errSeqOC[g, colIdx] <- as.numeric(scores[colIdx] > thr)
        
        
        
        if (updateMethod == "Simple") {
          alphaSeq[g, colIdx] <- alphaSeq[g, colIdx - 1] + gamma * (alpha - errSeqOC[g, colIdx])
          
        } else if (updateMethod == "Momentum") {
          w <- rev(momentumBW^(1:(colIdx - 1)))
          w <- w / sum(w)
          alphaSeq[g, colIdx] <- alphaSeq[g, colIdx - 1] + gamma * (alpha - sum(errSeqOC[g, colIdx] * w))
        }
        
        # Clipping
        alphaSeq[g, colIdx] <- min(upper_limit, max(lower_limit, alphaSeq[g, colIdx]))
      }
    }
  }
  
  
  cummean_cols <- apply(errSeqOC, 2, cummean)
  
  
  #  Aggregation with BOA
  aggalpha <- mixture(
    Y       = alphaSeq2, 
    experts = t(cummean_cols), 
    model   = "BOA", 
    loss.gradient = TRUE,
    loss.type = 'square'
  )
  
  print(aggalpha$weights)
  
  
  weights_df <- as.data.frame(aggalpha$weights)
  colnames(weights_df) <- paste0("Expert_", 1:K)  
  rownames(weights_df) <- paste0("T", 1:nrow(weights_df))  
  
  
  print(head(weights_df, 10))  
  
  
  for (t in T0:(myT - 1)) {
    colIdx <- t - T0 + 1
    
    historical <- scores[max(1, colIdx - lookback + 1):(colIdx - 1)]
    
    if (colIdx == 1) {
      aggregatedalpha[colIdx] <- alpha
      thr <- quantile(historical, 1 - aggregatedalpha[colIdx])
      
      err[colIdx] <- as.numeric(scores[colIdx] > thr)
      
    } else {
      w_t <- aggalpha$weights[colIdx - 1, ]
      aggregatedalpha[colIdx] <- sum(w_t * alphaSeq[, colIdx])
      
      thr <- quantile(historical, 1 - aggregatedalpha[colIdx])
      
      err[colIdx] <- as.numeric(scores[colIdx] > thr)
    }
  }
  
  
  # Restituzione dei risultati
  return(list(
    aggregatedalpha           = aggregatedalpha,
    err                       = err,
    errSeqOC                  = errSeqOC,
    garchForecastVec          = garchForecastVec,
    alphaSeq                  = alphaSeq,
    scores                    = scores,
    weights_df                = weights_df))
  
}

# === Application ===
getSymbols("AAPL", from = "2015-01-01", to = "2023-12-31")
## [1] "AAPL"
aapl_prices <- Cl(AAPL)
returns <- dailyReturn(aapl_prices)
returns <- na.omit(returns)


gammaGrid = c(0.03, 0.05, 0.07)
print(gammaGrid)
## [1] 0.03 0.05 0.07
results <- agaciVolatilityForecasting(
  returns   = returns, 
  alpha     = 0.1, 
  gammaGrid = gammaGrid, 
  lookback  = 1250, 
  garchP    = 1, 
  garchQ    = 1, 
  startUp   = 100,
  verbose   = TRUE, 
  updateMethod = "Simple", 
  momentumBW = 0.80
)
##                 X1         X2        X3
##    [1,] 0.33333333 0.33333333 0.3333333
##    [2,] 0.33333333 0.33333333 0.3333333
##    [3,] 0.33333333 0.33333333 0.3333333
##    [4,] 0.33333333 0.33333333 0.3333333
##    [5,] 0.33333333 0.33333333 0.3333333
##    [6,] 0.33333333 0.33333333 0.3333333
##    [7,] 0.33333333 0.33333333 0.3333333
##    [8,] 0.33333333 0.33333333 0.3333333
##    [9,] 0.33333333 0.33333333 0.3333333
##   [10,] 0.33333333 0.33333333 0.3333333
##   [11,] 0.33333333 0.33333333 0.3333333
##   [12,] 0.33333333 0.33333333 0.3333333
##   [13,] 0.33333333 0.33333333 0.3333333
##   [14,] 0.33333333 0.33333333 0.3333333
##   [15,] 0.33333333 0.33333333 0.3333333
##   [16,] 0.33333333 0.33333333 0.3333333
##   [17,] 0.33333333 0.33333333 0.3333333
##   [18,] 0.33333333 0.33333333 0.3333333
##   [19,] 0.33333333 0.33333333 0.3333333
##   [20,] 0.33333333 0.33333333 0.3333333
##   [21,] 0.33333333 0.33333333 0.3333333
##   [22,] 0.33333333 0.33333333 0.3333333
##   [23,] 0.33333333 0.33333333 0.3333333
##   [24,] 0.33333333 0.33333333 0.3333333
##   [25,] 0.33333333 0.33333333 0.3333333
##   [26,] 0.33333333 0.33333333 0.3333333
##   [27,] 0.33333333 0.33333333 0.3333333
##   [28,] 0.33333333 0.33333333 0.3333333
##   [29,] 0.33333333 0.33333333 0.3333333
##   [30,] 0.33333333 0.33333333 0.3333333
##   [31,] 0.33333333 0.33333333 0.3333333
##   [32,] 0.33333333 0.33333333 0.3333333
##   [33,] 0.33333333 0.33333333 0.3333333
##   [34,] 0.33333333 0.33333333 0.3333333
##   [35,] 0.33333333 0.33333333 0.3333333
##   [36,] 0.33333333 0.33333333 0.3333333
##   [37,] 0.33333333 0.33333333 0.3333333
##   [38,] 0.33333333 0.33333333 0.3333333
##   [39,] 0.33333333 0.33333333 0.3333333
##   [40,] 0.33333333 0.33333333 0.3333333
##   [41,] 0.33333333 0.33333333 0.3333333
##   [42,] 0.33333333 0.33333333 0.3333333
##   [43,] 0.33333333 0.33333333 0.3333333
##   [44,] 0.33333333 0.33333333 0.3333333
##   [45,] 0.33333333 0.33333333 0.3333333
##   [46,] 0.46993293 0.36928608 0.1607810
##   [47,] 0.46993293 0.36928608 0.1607810
##   [48,] 0.46993293 0.36928608 0.1607810
##   [49,] 0.16873597 0.35744152 0.4738225
##   [50,] 0.16873597 0.35744152 0.4738225
##   [51,] 0.15056007 0.31893863 0.5305013
##   [52,] 0.15056007 0.31893863 0.5305013
##   [53,] 0.15056007 0.31893863 0.5305013
##   [54,] 0.15056007 0.31893863 0.5305013
##   [55,] 0.15056007 0.31893863 0.5305013
##   [56,] 0.15056007 0.31893863 0.5305013
##   [57,] 0.15056007 0.31893863 0.5305013
##   [58,] 0.15056007 0.31893863 0.5305013
##   [59,] 0.15463120 0.32756271 0.5178061
##   [60,] 0.15463120 0.32756271 0.5178061
##   [61,] 0.11562524 0.24493451 0.6394403
##   [62,] 0.11562524 0.24493451 0.6394403
##   [63,] 0.11562524 0.24493451 0.6394403
##   [64,] 0.11562524 0.24493451 0.6394403
##   [65,] 0.11562524 0.24493451 0.6394403
##   [66,] 0.11562524 0.24493451 0.6394403
##   [67,] 0.11562524 0.24493451 0.6394403
##   [68,] 0.11562524 0.24493451 0.6394403
##   [69,] 0.11562524 0.24493451 0.6394403
##   [70,] 0.11562524 0.24493451 0.6394403
##   [71,] 0.11562524 0.24493451 0.6394403
##   [72,] 0.11562524 0.24493451 0.6394403
##   [73,] 0.11562524 0.24493451 0.6394403
##   [74,] 0.11562524 0.24493451 0.6394403
##   [75,] 0.11562524 0.24493451 0.6394403
##   [76,] 0.17518385 0.25722522 0.5675909
##   [77,] 0.17518385 0.25722522 0.5675909
##   [78,] 0.17518385 0.25722522 0.5675909
##   [79,] 0.17518385 0.25722522 0.5675909
##   [80,] 0.17518385 0.25722522 0.5675909
##   [81,] 0.17518385 0.25722522 0.5675909
##   [82,] 0.17518385 0.25722522 0.5675909
##   [83,] 0.17518385 0.25722522 0.5675909
##   [84,] 0.17518385 0.25722522 0.5675909
##   [85,] 0.17518385 0.25722522 0.5675909
##   [86,] 0.17518385 0.25722522 0.5675909
##   [87,] 0.17518385 0.25722522 0.5675909
##   [88,] 0.17518385 0.25722522 0.5675909
##   [89,] 0.17518385 0.25722522 0.5675909
##   [90,] 0.17518385 0.25722522 0.5675909
##   [91,] 0.17518385 0.25722522 0.5675909
##   [92,] 0.17518385 0.25722522 0.5675909
##   [93,] 0.17518385 0.25722522 0.5675909
##   [94,] 0.17518385 0.25722522 0.5675909
##   [95,] 0.17518385 0.25722522 0.5675909
##   [96,] 0.17518385 0.25722522 0.5675909
##   [97,] 0.17518385 0.25722522 0.5675909
##   [98,] 0.17518385 0.25722522 0.5675909
##   [99,] 0.17518385 0.25722522 0.5675909
##  [100,] 0.17518385 0.25722522 0.5675909
##  [101,] 0.17518385 0.25722522 0.5675909
##  [102,] 0.17518385 0.25722522 0.5675909
##  [103,] 0.17518385 0.25722522 0.5675909
##  [104,] 0.17518385 0.25722522 0.5675909
##  [105,] 0.17971443 0.26622355 0.5540620
##  [106,] 0.17971443 0.26622355 0.5540620
##  [107,] 0.17971443 0.26622355 0.5540620
##  [108,] 0.17971443 0.26622355 0.5540620
##  [109,] 0.17971443 0.26622355 0.5540620
##  [110,] 0.17971443 0.26622355 0.5540620
##  [111,] 0.17971443 0.26622355 0.5540620
##  [112,] 0.17971443 0.26622355 0.5540620
##  [113,] 0.17971443 0.26622355 0.5540620
##  [114,] 0.17971443 0.26622355 0.5540620
##  [115,] 0.17971443 0.26622355 0.5540620
##  [116,] 0.17971443 0.26622355 0.5540620
##  [117,] 0.17971443 0.26622355 0.5540620
##  [118,] 0.17971443 0.26622355 0.5540620
##  [119,] 0.17971443 0.26622355 0.5540620
##  [120,] 0.17971443 0.26622355 0.5540620
##  [121,] 0.17971443 0.26622355 0.5540620
##  [122,] 0.17971443 0.26622355 0.5540620
##  [123,] 0.17971443 0.26622355 0.5540620
##  [124,] 0.17971443 0.26622355 0.5540620
##  [125,] 0.17971443 0.26622355 0.5540620
##  [126,] 0.17971443 0.26622355 0.5540620
##  [127,] 0.17971443 0.26622355 0.5540620
##  [128,] 0.17971443 0.26622355 0.5540620
##  [129,] 0.17971443 0.26622355 0.5540620
##  [130,] 0.17971443 0.26622355 0.5540620
##  [131,] 0.17971443 0.26622355 0.5540620
##  [132,] 0.17971443 0.26622355 0.5540620
##  [133,] 0.17971443 0.26622355 0.5540620
##  [134,] 0.17971443 0.26622355 0.5540620
##  [135,] 0.17971443 0.26622355 0.5540620
##  [136,] 0.17971443 0.26622355 0.5540620
##  [137,] 0.17971443 0.26622355 0.5540620
##  [138,] 0.17971443 0.26622355 0.5540620
##  [139,] 0.17971443 0.26622355 0.5540620
##  [140,] 0.17971443 0.26622355 0.5540620
##  [141,] 0.17971443 0.26622355 0.5540620
##  [142,] 0.17971443 0.26622355 0.5540620
##  [143,] 0.17971443 0.26622355 0.5540620
##  [144,] 0.17971443 0.26622355 0.5540620
##  [145,] 0.17971443 0.26622355 0.5540620
##  [146,] 0.17971443 0.26622355 0.5540620
##  [147,] 0.17971443 0.26622355 0.5540620
##  [148,] 0.17971443 0.26622355 0.5540620
##  [149,] 0.17971443 0.26622355 0.5540620
##  [150,] 0.17971443 0.26622355 0.5540620
##  [151,] 0.17971443 0.26622355 0.5540620
##  [152,] 0.17971443 0.26622355 0.5540620
##  [153,] 0.17971443 0.26622355 0.5540620
##  [154,] 0.17971443 0.26622355 0.5540620
##  [155,] 0.17971443 0.26622355 0.5540620
##  [156,] 0.17971443 0.26622355 0.5540620
##  [157,] 0.17971443 0.26622355 0.5540620
##  [158,] 0.17971443 0.26622355 0.5540620
##  [159,] 0.17971443 0.26622355 0.5540620
##  [160,] 0.17971443 0.26622355 0.5540620
##  [161,] 0.25372870 0.27064003 0.4756313
##  [162,] 0.25372870 0.27064003 0.4756313
##  [163,] 0.25372870 0.27064003 0.4756313
##  [164,] 0.25372870 0.27064003 0.4756313
##  [165,] 0.25372870 0.27064003 0.4756313
##  [166,] 0.25372870 0.27064003 0.4756313
##  [167,] 0.25372870 0.27064003 0.4756313
##  [168,] 0.25372870 0.27064003 0.4756313
##  [169,] 0.25372870 0.27064003 0.4756313
##  [170,] 0.25372870 0.27064003 0.4756313
##  [171,] 0.25372870 0.27064003 0.4756313
##  [172,] 0.25372870 0.27064003 0.4756313
##  [173,] 0.25372870 0.27064003 0.4756313
##  [174,] 0.25372870 0.27064003 0.4756313
##  [175,] 0.25372870 0.27064003 0.4756313
##  [176,] 0.25372870 0.27064003 0.4756313
##  [177,] 0.25372870 0.27064003 0.4756313
##  [178,] 0.25372870 0.27064003 0.4756313
##  [179,] 0.25372870 0.27064003 0.4756313
##  [180,] 0.25372870 0.27064003 0.4756313
##  [181,] 0.25372870 0.27064003 0.4756313
##  [182,] 0.25372870 0.27064003 0.4756313
##  [183,] 0.25372870 0.27064003 0.4756313
##  [184,] 0.25372870 0.27064003 0.4756313
##  [185,] 0.25372870 0.27064003 0.4756313
##  [186,] 0.25372870 0.27064003 0.4756313
##  [187,] 0.25372870 0.27064003 0.4756313
##  [188,] 0.25372870 0.27064003 0.4756313
##  [189,] 0.25372870 0.27064003 0.4756313
##  [190,] 0.25372870 0.27064003 0.4756313
##  [191,] 0.25372870 0.27064003 0.4756313
##  [192,] 0.25372870 0.27064003 0.4756313
##  [193,] 0.25372870 0.27064003 0.4756313
##  [194,] 0.25372870 0.27064003 0.4756313
##  [195,] 0.25372870 0.27064003 0.4756313
##  [196,] 0.25372870 0.27064003 0.4756313
##  [197,] 0.25372870 0.27064003 0.4756313
##  [198,] 0.25372870 0.27064003 0.4756313
##  [199,] 0.25372870 0.27064003 0.4756313
##  [200,] 0.25372870 0.27064003 0.4756313
##  [201,] 0.25372870 0.27064003 0.4756313
##  [202,] 0.25372870 0.27064003 0.4756313
##  [203,] 0.25372870 0.27064003 0.4756313
##  [204,] 0.25372870 0.27064003 0.4756313
##  [205,] 0.25372870 0.27064003 0.4756313
##  [206,] 0.25372870 0.27064003 0.4756313
##  [207,] 0.25372870 0.27064003 0.4756313
##  [208,] 0.25372870 0.27064003 0.4756313
##  [209,] 0.25372870 0.27064003 0.4756313
##  [210,] 0.25372870 0.27064003 0.4756313
##  [211,] 0.25372870 0.27064003 0.4756313
##  [212,] 0.25372870 0.27064003 0.4756313
##  [213,] 0.25372870 0.27064003 0.4756313
##  [214,] 0.25372870 0.27064003 0.4756313
##  [215,] 0.25372870 0.27064003 0.4756313
##  [216,] 0.25372870 0.27064003 0.4756313
##  [217,] 0.25372870 0.27064003 0.4756313
##  [218,] 0.25372870 0.27064003 0.4756313
##  [219,] 0.25372870 0.27064003 0.4756313
##  [220,] 0.25372870 0.27064003 0.4756313
##  [221,] 0.25717940 0.27615067 0.4666699
##  [222,] 0.25717940 0.27615067 0.4666699
##  [223,] 0.25717940 0.27615067 0.4666699
##  [224,] 0.25717940 0.27615067 0.4666699
##  [225,] 0.25717940 0.27615067 0.4666699
##  [226,] 0.25717940 0.27615067 0.4666699
##  [227,] 0.25717940 0.27615067 0.4666699
##  [228,] 0.25717940 0.27615067 0.4666699
##  [229,] 0.25717940 0.27615067 0.4666699
##  [230,] 0.25717940 0.27615067 0.4666699
##  [231,] 0.25717940 0.27615067 0.4666699
##  [232,] 0.25717940 0.27615067 0.4666699
##  [233,] 0.25717940 0.27615067 0.4666699
##  [234,] 0.25717940 0.27615067 0.4666699
##  [235,] 0.25717940 0.27615067 0.4666699
##  [236,] 0.25717940 0.27615067 0.4666699
##  [237,] 0.25717940 0.27615067 0.4666699
##  [238,] 0.25717940 0.27615067 0.4666699
##  [239,] 0.25717940 0.27615067 0.4666699
##  [240,] 0.25717940 0.27615067 0.4666699
##  [241,] 0.25717940 0.27615067 0.4666699
##  [242,] 0.16091682 0.15194304 0.6871401
##  [243,] 0.16091682 0.15194304 0.6871401
##  [244,] 0.16091682 0.15194304 0.6871401
##  [245,] 0.16091682 0.15194304 0.6871401
##  [246,] 0.16091682 0.15194304 0.6871401
##  [247,] 0.16091682 0.15194304 0.6871401
##  [248,] 0.16091682 0.15194304 0.6871401
##  [249,] 0.16091682 0.15194304 0.6871401
##  [250,] 0.16091682 0.15194304 0.6871401
##  [251,] 0.16091682 0.15194304 0.6871401
##  [252,] 0.16091682 0.15194304 0.6871401
##  [253,] 0.16091682 0.15194304 0.6871401
##  [254,] 0.16091682 0.15194304 0.6871401
##  [255,] 0.16091682 0.15194304 0.6871401
##  [256,] 0.16091682 0.15194304 0.6871401
##  [257,] 0.16713305 0.16028513 0.6725818
##  [258,] 0.16713305 0.16028513 0.6725818
##  [259,] 0.16713305 0.16028513 0.6725818
##  [260,] 0.16713305 0.16028513 0.6725818
##  [261,] 0.16713305 0.16028513 0.6725818
##  [262,] 0.16713305 0.16028513 0.6725818
##  [263,] 0.16713305 0.16028513 0.6725818
##  [264,] 0.16713305 0.16028513 0.6725818
##  [265,] 0.16713305 0.16028513 0.6725818
##  [266,] 0.16713305 0.16028513 0.6725818
##  [267,] 0.16713305 0.16028513 0.6725818
##  [268,] 0.16713305 0.16028513 0.6725818
##  [269,] 0.16713305 0.16028513 0.6725818
##  [270,] 0.16713305 0.16028513 0.6725818
##  [271,] 0.16713305 0.16028513 0.6725818
##  [272,] 0.16713305 0.16028513 0.6725818
##  [273,] 0.16713305 0.16028513 0.6725818
##  [274,] 0.16713305 0.16028513 0.6725818
##  [275,] 0.16713305 0.16028513 0.6725818
##  [276,] 0.16713305 0.16028513 0.6725818
##  [277,] 0.16713305 0.16028513 0.6725818
##  [278,] 0.16713305 0.16028513 0.6725818
##  [279,] 0.16713305 0.16028513 0.6725818
##  [280,] 0.16713305 0.16028513 0.6725818
##  [281,] 0.16713305 0.16028513 0.6725818
##  [282,] 0.16713305 0.16028513 0.6725818
##  [283,] 0.16713305 0.16028513 0.6725818
##  [284,] 0.16713305 0.16028513 0.6725818
##  [285,] 0.16713305 0.16028513 0.6725818
##  [286,] 0.16713305 0.16028513 0.6725818
##  [287,] 0.16713305 0.16028513 0.6725818
##  [288,] 0.16713305 0.16028513 0.6725818
##  [289,] 0.16713305 0.16028513 0.6725818
##  [290,] 0.16713305 0.16028513 0.6725818
##  [291,] 0.16713305 0.16028513 0.6725818
##  [292,] 0.16713305 0.16028513 0.6725818
##  [293,] 0.16713305 0.16028513 0.6725818
##  [294,] 0.16713305 0.16028513 0.6725818
##  [295,] 0.16713305 0.16028513 0.6725818
##  [296,] 0.16713305 0.16028513 0.6725818
##  [297,] 0.16713305 0.16028513 0.6725818
##  [298,] 0.16713305 0.16028513 0.6725818
##  [299,] 0.16713305 0.16028513 0.6725818
##  [300,] 0.13184392 0.10710749 0.7610486
##  [301,] 0.13184392 0.10710749 0.7610486
##  [302,] 0.13184392 0.10710749 0.7610486
##  [303,] 0.13184392 0.10710749 0.7610486
##  [304,] 0.13184392 0.10710749 0.7610486
##  [305,] 0.13184392 0.10710749 0.7610486
##  [306,] 0.13184392 0.10710749 0.7610486
##  [307,] 0.13184392 0.10710749 0.7610486
##  [308,] 0.13184392 0.10710749 0.7610486
##  [309,] 0.13184392 0.10710749 0.7610486
##  [310,] 0.13184392 0.10710749 0.7610486
##  [311,] 0.13184392 0.10710749 0.7610486
##  [312,] 0.13184392 0.10710749 0.7610486
##  [313,] 0.13184392 0.10710749 0.7610486
##  [314,] 0.13184392 0.10710749 0.7610486
##  [315,] 0.13820419 0.11411571 0.7476801
##  [316,] 0.13820419 0.11411571 0.7476801
##  [317,] 0.13820419 0.11411571 0.7476801
##  [318,] 0.13820419 0.11411571 0.7476801
##  [319,] 0.13820419 0.11411571 0.7476801
##  [320,] 0.13820419 0.11411571 0.7476801
##  [321,] 0.13820419 0.11411571 0.7476801
##  [322,] 0.13820419 0.11411571 0.7476801
##  [323,] 0.13820419 0.11411571 0.7476801
##  [324,] 0.13820419 0.11411571 0.7476801
##  [325,] 0.13820419 0.11411571 0.7476801
##  [326,] 0.13820419 0.11411571 0.7476801
##  [327,] 0.13820419 0.11411571 0.7476801
##  [328,] 0.13820419 0.11411571 0.7476801
##  [329,] 0.13820419 0.11411571 0.7476801
##  [330,] 0.13820419 0.11411571 0.7476801
##  [331,] 0.13820419 0.11411571 0.7476801
##  [332,] 0.22980711 0.11995221 0.6502407
##  [333,] 0.22980711 0.11995221 0.6502407
##  [334,] 0.22980711 0.11995221 0.6502407
##  [335,] 0.22980711 0.11995221 0.6502407
##  [336,] 0.22980711 0.11995221 0.6502407
##  [337,] 0.22980711 0.11995221 0.6502407
##  [338,] 0.22980711 0.11995221 0.6502407
##  [339,] 0.22980711 0.11995221 0.6502407
##  [340,] 0.22980711 0.11995221 0.6502407
##  [341,] 0.22980711 0.11995221 0.6502407
##  [342,] 0.22980711 0.11995221 0.6502407
##  [343,] 0.22980711 0.11995221 0.6502407
##  [344,] 0.22980711 0.11995221 0.6502407
##  [345,] 0.22980711 0.11995221 0.6502407
##  [346,] 0.22980711 0.11995221 0.6502407
##  [347,] 0.22980711 0.11995221 0.6502407
##  [348,] 0.22980711 0.11995221 0.6502407
##  [349,] 0.22980711 0.11995221 0.6502407
##  [350,] 0.22980711 0.11995221 0.6502407
##  [351,] 0.17520089 0.13305784 0.6917413
##  [352,] 0.17520089 0.13305784 0.6917413
##  [353,] 0.17520089 0.13305784 0.6917413
##  [354,] 0.17520089 0.13305784 0.6917413
##  [355,] 0.17520089 0.13305784 0.6917413
##  [356,] 0.17520089 0.13305784 0.6917413
##  [357,] 0.17520089 0.13305784 0.6917413
##  [358,] 0.17520089 0.13305784 0.6917413
##  [359,] 0.17520089 0.13305784 0.6917413
##  [360,] 0.17520089 0.13305784 0.6917413
##  [361,] 0.17520089 0.13305784 0.6917413
##  [362,] 0.17520089 0.13305784 0.6917413
##  [363,] 0.18486829 0.12539473 0.6897370
##  [364,] 0.18486829 0.12539473 0.6897370
##  [365,] 0.18486829 0.12539473 0.6897370
##  [366,] 0.18486829 0.12539473 0.6897370
##  [367,] 0.18486829 0.12539473 0.6897370
##  [368,] 0.18486829 0.12539473 0.6897370
##  [369,] 0.18486829 0.12539473 0.6897370
##  [370,] 0.18486829 0.12539473 0.6897370
##  [371,] 0.18486829 0.12539473 0.6897370
##  [372,] 0.18486829 0.12539473 0.6897370
##  [373,] 0.18486829 0.12539473 0.6897370
##  [374,] 0.18486829 0.12539473 0.6897370
##  [375,] 0.18486829 0.12539473 0.6897370
##  [376,] 0.18486829 0.12539473 0.6897370
##  [377,] 0.18486829 0.12539473 0.6897370
##  [378,] 0.18486829 0.12539473 0.6897370
##  [379,] 0.18486829 0.12539473 0.6897370
##  [380,] 0.18486829 0.12539473 0.6897370
##  [381,] 0.18486829 0.12539473 0.6897370
##  [382,] 0.18486829 0.12539473 0.6897370
##  [383,] 0.18486829 0.12539473 0.6897370
##  [384,] 0.18486829 0.12539473 0.6897370
##  [385,] 0.18486829 0.12539473 0.6897370
##  [386,] 0.18486829 0.12539473 0.6897370
##  [387,] 0.18486829 0.12539473 0.6897370
##  [388,] 0.18486829 0.12539473 0.6897370
##  [389,] 0.18486829 0.12539473 0.6897370
##  [390,] 0.18486829 0.12539473 0.6897370
##  [391,] 0.18486829 0.12539473 0.6897370
##  [392,] 0.18486829 0.12539473 0.6897370
##  [393,] 0.18486829 0.12539473 0.6897370
##  [394,] 0.18486829 0.12539473 0.6897370
##  [395,] 0.18486829 0.12539473 0.6897370
##  [396,] 0.18486829 0.12539473 0.6897370
##  [397,] 0.18486829 0.12539473 0.6897370
##  [398,] 0.18486829 0.12539473 0.6897370
##  [399,] 0.18486829 0.12539473 0.6897370
##  [400,] 0.18486829 0.12539473 0.6897370
##  [401,] 0.18486829 0.12539473 0.6897370
##  [402,] 0.18486829 0.12539473 0.6897370
##  [403,] 0.18486829 0.12539473 0.6897370
##  [404,] 0.18486829 0.12539473 0.6897370
##  [405,] 0.18486829 0.12539473 0.6897370
##  [406,] 0.18486829 0.12539473 0.6897370
##  [407,] 0.18486829 0.12539473 0.6897370
##  [408,] 0.18486829 0.12539473 0.6897370
##  [409,] 0.18486829 0.12539473 0.6897370
##  [410,] 0.18486829 0.12539473 0.6897370
##  [411,] 0.18486829 0.12539473 0.6897370
##  [412,] 0.18486829 0.12539473 0.6897370
##  [413,] 0.18486829 0.12539473 0.6897370
##  [414,] 0.18486829 0.12539473 0.6897370
##  [415,] 0.18486829 0.12539473 0.6897370
##  [416,] 0.18486829 0.12539473 0.6897370
##  [417,] 0.18486829 0.12539473 0.6897370
##  [418,] 0.18486829 0.12539473 0.6897370
##  [419,] 0.18486829 0.12539473 0.6897370
##  [420,] 0.18486829 0.12539473 0.6897370
##  [421,] 0.18486829 0.12539473 0.6897370
##  [422,] 0.18486829 0.12539473 0.6897370
##  [423,] 0.18486829 0.12539473 0.6897370
##  [424,] 0.18486829 0.12539473 0.6897370
##  [425,] 0.18486829 0.12539473 0.6897370
##  [426,] 0.18486829 0.12539473 0.6897370
##  [427,] 0.18486829 0.12539473 0.6897370
##  [428,] 0.18486829 0.12539473 0.6897370
##  [429,] 0.18486829 0.12539473 0.6897370
##  [430,] 0.18486829 0.12539473 0.6897370
##  [431,] 0.18486829 0.12539473 0.6897370
##  [432,] 0.18486829 0.12539473 0.6897370
##  [433,] 0.18486829 0.12539473 0.6897370
##  [434,] 0.18486829 0.12539473 0.6897370
##  [435,] 0.18486829 0.12539473 0.6897370
##  [436,] 0.18486829 0.12539473 0.6897370
##  [437,] 0.18486829 0.12539473 0.6897370
##  [438,] 0.18486829 0.12539473 0.6897370
##  [439,] 0.18486829 0.12539473 0.6897370
##  [440,] 0.18486829 0.12539473 0.6897370
##  [441,] 0.18486829 0.12539473 0.6897370
##  [442,] 0.18486829 0.12539473 0.6897370
##  [443,] 0.18486829 0.12539473 0.6897370
##  [444,] 0.18486829 0.12539473 0.6897370
##  [445,] 0.18486829 0.12539473 0.6897370
##  [446,] 0.18486829 0.12539473 0.6897370
##  [447,] 0.18486829 0.12539473 0.6897370
##  [448,] 0.18486829 0.12539473 0.6897370
##  [449,] 0.18486829 0.12539473 0.6897370
##  [450,] 0.18486829 0.12539473 0.6897370
##  [451,] 0.18486829 0.12539473 0.6897370
##  [452,] 0.18486829 0.12539473 0.6897370
##  [453,] 0.18486829 0.12539473 0.6897370
##  [454,] 0.18486829 0.12539473 0.6897370
##  [455,] 0.18486829 0.12539473 0.6897370
##  [456,] 0.18486829 0.12539473 0.6897370
##  [457,] 0.18486829 0.12539473 0.6897370
##  [458,] 0.18486829 0.12539473 0.6897370
##  [459,] 0.18486829 0.12539473 0.6897370
##  [460,] 0.18486829 0.12539473 0.6897370
##  [461,] 0.18486829 0.12539473 0.6897370
##  [462,] 0.18486829 0.12539473 0.6897370
##  [463,] 0.18486829 0.12539473 0.6897370
##  [464,] 0.18486829 0.12539473 0.6897370
##  [465,] 0.18486829 0.12539473 0.6897370
##  [466,] 0.18486829 0.12539473 0.6897370
##  [467,] 0.18486829 0.12539473 0.6897370
##  [468,] 0.18486829 0.12539473 0.6897370
##  [469,] 0.18486829 0.12539473 0.6897370
##  [470,] 0.18486829 0.12539473 0.6897370
##  [471,] 0.18486829 0.12539473 0.6897370
##  [472,] 0.13730668 0.12649923 0.7361941
##  [473,] 0.13730668 0.12649923 0.7361941
##  [474,] 0.13730668 0.12649923 0.7361941
##  [475,] 0.13730668 0.12649923 0.7361941
##  [476,] 0.13730668 0.12649923 0.7361941
##  [477,] 0.13730668 0.12649923 0.7361941
##  [478,] 0.13730668 0.12649923 0.7361941
##  [479,] 0.13730668 0.12649923 0.7361941
##  [480,] 0.13730668 0.12649923 0.7361941
##  [481,] 0.13730668 0.12649923 0.7361941
##  [482,] 0.13730668 0.12649923 0.7361941
##  [483,] 0.13730668 0.12649923 0.7361941
##  [484,] 0.13730668 0.12649923 0.7361941
##  [485,] 0.13730668 0.12649923 0.7361941
##  [486,] 0.13730668 0.12649923 0.7361941
##  [487,] 0.13730668 0.12649923 0.7361941
##  [488,] 0.13730668 0.12649923 0.7361941
##  [489,] 0.13730668 0.12649923 0.7361941
##  [490,] 0.13730668 0.12649923 0.7361941
##  [491,] 0.13730668 0.12649923 0.7361941
##  [492,] 0.13730668 0.12649923 0.7361941
##  [493,] 0.11451545 0.09093155 0.7945530
##  [494,] 0.11451545 0.09093155 0.7945530
##  [495,] 0.11451545 0.09093155 0.7945530
##  [496,] 0.11451545 0.09093155 0.7945530
##  [497,] 0.11451545 0.09093155 0.7945530
##  [498,] 0.11451545 0.09093155 0.7945530
##  [499,] 0.11451545 0.09093155 0.7945530
##  [500,] 0.11451545 0.09093155 0.7945530
##  [501,] 0.11451545 0.09093155 0.7945530
##  [502,] 0.11451545 0.09093155 0.7945530
##  [503,] 0.11451545 0.09093155 0.7945530
##  [504,] 0.11451545 0.09093155 0.7945530
##  [505,] 0.09022057 0.09966048 0.8101190
##  [506,] 0.09022057 0.09966048 0.8101190
##  [507,] 0.09022057 0.09966048 0.8101190
##  [508,] 0.09022057 0.09966048 0.8101190
##  [509,] 0.09022057 0.09966048 0.8101190
##  [510,] 0.09022057 0.09966048 0.8101190
##  [511,] 0.09022057 0.09966048 0.8101190
##  [512,] 0.09022057 0.09966048 0.8101190
##  [513,] 0.09022057 0.09966048 0.8101190
##  [514,] 0.09022057 0.09966048 0.8101190
##  [515,] 0.09022057 0.09966048 0.8101190
##  [516,] 0.09022057 0.09966048 0.8101190
##  [517,] 0.09022057 0.09966048 0.8101190
##  [518,] 0.09022057 0.09966048 0.8101190
##  [519,] 0.09022057 0.09966048 0.8101190
##  [520,] 0.09022057 0.09966048 0.8101190
##  [521,] 0.09022057 0.09966048 0.8101190
##  [522,] 0.09022057 0.09966048 0.8101190
##  [523,] 0.09022057 0.09966048 0.8101190
##  [524,] 0.09022057 0.09966048 0.8101190
##  [525,] 0.09022057 0.09966048 0.8101190
##  [526,] 0.09022057 0.09966048 0.8101190
##  [527,] 0.09022057 0.09966048 0.8101190
##  [528,] 0.09022057 0.09966048 0.8101190
##  [529,] 0.09022057 0.09966048 0.8101190
##  [530,] 0.09022057 0.09966048 0.8101190
##  [531,] 0.09022057 0.09966048 0.8101190
##  [532,] 0.09022057 0.09966048 0.8101190
##  [533,] 0.09022057 0.09966048 0.8101190
##  [534,] 0.09022057 0.09966048 0.8101190
##  [535,] 0.09022057 0.09966048 0.8101190
##  [536,] 0.09022057 0.09966048 0.8101190
##  [537,] 0.09022057 0.09966048 0.8101190
##  [538,] 0.09022057 0.09966048 0.8101190
##  [539,] 0.09022057 0.09966048 0.8101190
##  [540,] 0.09022057 0.09966048 0.8101190
##  [541,] 0.09022057 0.09966048 0.8101190
##  [542,] 0.09022057 0.09966048 0.8101190
##  [543,] 0.10164570 0.11176221 0.7865921
##  [544,] 0.10164570 0.11176221 0.7865921
##  [545,] 0.10164570 0.11176221 0.7865921
##  [546,] 0.10164570 0.11176221 0.7865921
##  [547,] 0.10164570 0.11176221 0.7865921
##  [548,] 0.10164570 0.11176221 0.7865921
##  [549,] 0.10164570 0.11176221 0.7865921
##  [550,] 0.10164570 0.11176221 0.7865921
##  [551,] 0.10164570 0.11176221 0.7865921
##  [552,] 0.10164570 0.11176221 0.7865921
##  [553,] 0.10164570 0.11176221 0.7865921
##  [554,] 0.10164570 0.11176221 0.7865921
##  [555,] 0.10164570 0.11176221 0.7865921
##  [556,] 0.10164570 0.11176221 0.7865921
##  [557,] 0.10164570 0.11176221 0.7865921
##  [558,] 0.10164570 0.11176221 0.7865921
##  [559,] 0.10164570 0.11176221 0.7865921
##  [560,] 0.10164570 0.11176221 0.7865921
##  [561,] 0.10164570 0.11176221 0.7865921
##  [562,] 0.10164570 0.11176221 0.7865921
##  [563,] 0.10164570 0.11176221 0.7865921
##  [564,] 0.10164570 0.11176221 0.7865921
##  [565,] 0.10164570 0.11176221 0.7865921
##  [566,] 0.10164570 0.11176221 0.7865921
##  [567,] 0.10164570 0.11176221 0.7865921
##  [568,] 0.10164570 0.11176221 0.7865921
##  [569,] 0.10164570 0.11176221 0.7865921
##  [570,] 0.10164570 0.11176221 0.7865921
##  [571,] 0.10164570 0.11176221 0.7865921
##  [572,] 0.10164570 0.11176221 0.7865921
##  [573,] 0.10164570 0.11176221 0.7865921
##  [574,] 0.10164570 0.11176221 0.7865921
##  [575,] 0.10164570 0.11176221 0.7865921
##  [576,] 0.10164570 0.11176221 0.7865921
##  [577,] 0.10164570 0.11176221 0.7865921
##  [578,] 0.11517511 0.12529240 0.7595325
##  [579,] 0.11517511 0.12529240 0.7595325
##  [580,] 0.11517511 0.12529240 0.7595325
##  [581,] 0.11517511 0.12529240 0.7595325
##  [582,] 0.11921692 0.13212773 0.7486554
##  [583,] 0.11921692 0.13212773 0.7486554
##  [584,] 0.11921692 0.13212773 0.7486554
##  [585,] 0.11921692 0.13212773 0.7486554
##  [586,] 0.11921692 0.13212773 0.7486554
##  [587,] 0.11921692 0.13212773 0.7486554
##  [588,] 0.11921692 0.13212773 0.7486554
##  [589,] 0.11921692 0.13212773 0.7486554
##  [590,] 0.10364430 0.09877382 0.7975819
##  [591,] 0.10364430 0.09877382 0.7975819
##  [592,] 0.10364430 0.09877382 0.7975819
##  [593,] 0.10364430 0.09877382 0.7975819
##  [594,] 0.10364430 0.09877382 0.7975819
##  [595,] 0.10364430 0.09877382 0.7975819
##  [596,] 0.10364430 0.09877382 0.7975819
##  [597,] 0.10364430 0.09877382 0.7975819
##  [598,] 0.10364430 0.09877382 0.7975819
##  [599,] 0.10364430 0.09877382 0.7975819
##  [600,] 0.08603384 0.09475399 0.8192122
##  [601,] 0.08603384 0.09475399 0.8192122
##  [602,] 0.08603384 0.09475399 0.8192122
##  [603,] 0.08603384 0.09475399 0.8192122
##  [604,] 0.08603384 0.09475399 0.8192122
##  [605,] 0.07236896 0.09050363 0.8371274
##  [606,] 0.07236896 0.09050363 0.8371274
##  [607,] 0.07236896 0.09050363 0.8371274
##  [608,] 0.07236896 0.09050363 0.8371274
##  [609,] 0.07236896 0.09050363 0.8371274
##  [610,] 0.06243083 0.06809132 0.8694778
##  [611,] 0.06243083 0.06809132 0.8694778
##  [612,] 0.06243083 0.06809132 0.8694778
##  [613,] 0.06243083 0.06809132 0.8694778
##  [614,] 0.06243083 0.06809132 0.8694778
##  [615,] 0.06243083 0.06809132 0.8694778
##  [616,] 0.06243083 0.06809132 0.8694778
##  [617,] 0.06243083 0.06809132 0.8694778
##  [618,] 0.06243083 0.06809132 0.8694778
##  [619,] 0.06243083 0.06809132 0.8694778
##  [620,] 0.06243083 0.06809132 0.8694778
##  [621,] 0.06243083 0.06809132 0.8694778
##  [622,] 0.06243083 0.06809132 0.8694778
##  [623,] 0.06243083 0.06809132 0.8694778
##  [624,] 0.06243083 0.06809132 0.8694778
##  [625,] 0.06243083 0.06809132 0.8694778
##  [626,] 0.06243083 0.06809132 0.8694778
##  [627,] 0.06243083 0.06809132 0.8694778
##  [628,] 0.06243083 0.06809132 0.8694778
##  [629,] 0.06243083 0.06809132 0.8694778
##  [630,] 0.06243083 0.06809132 0.8694778
##  [631,] 0.06243083 0.06809132 0.8694778
##  [632,] 0.06243083 0.06809132 0.8694778
##  [633,] 0.06243083 0.06809132 0.8694778
##  [634,] 0.06243083 0.06809132 0.8694778
##  [635,] 0.06243083 0.06809132 0.8694778
##  [636,] 0.06243083 0.06809132 0.8694778
##  [637,] 0.06243083 0.06809132 0.8694778
##  [638,] 0.06243083 0.06809132 0.8694778
##  [639,] 0.06243083 0.06809132 0.8694778
##  [640,] 0.06243083 0.06809132 0.8694778
##  [641,] 0.06243083 0.06809132 0.8694778
##  [642,] 0.06243083 0.06809132 0.8694778
##  [643,] 0.06243083 0.06809132 0.8694778
##  [644,] 0.06243083 0.06809132 0.8694778
##  [645,] 0.06243083 0.06809132 0.8694778
##  [646,] 0.06243083 0.06809132 0.8694778
##  [647,] 0.06243083 0.06809132 0.8694778
##  [648,] 0.06243083 0.06809132 0.8694778
##  [649,] 0.06243083 0.06809132 0.8694778
##  [650,] 0.06243083 0.06809132 0.8694778
##  [651,] 0.07078401 0.07511903 0.8540970
##  [652,] 0.07078401 0.07511903 0.8540970
##  [653,] 0.07078401 0.07511903 0.8540970
##  [654,] 0.07078401 0.07511903 0.8540970
##  [655,] 0.07078401 0.07511903 0.8540970
##  [656,] 0.07078401 0.07511903 0.8540970
##  [657,] 0.07078401 0.07511903 0.8540970
##  [658,] 0.07078401 0.07511903 0.8540970
##  [659,] 0.07078401 0.07511903 0.8540970
##  [660,] 0.07078401 0.07511903 0.8540970
##  [661,] 0.07078401 0.07511903 0.8540970
##  [662,] 0.07078401 0.07511903 0.8540970
##  [663,] 0.07078401 0.07511903 0.8540970
##  [664,] 0.07078401 0.07511903 0.8540970
##  [665,] 0.07078401 0.07511903 0.8540970
##  [666,] 0.07078401 0.07511903 0.8540970
##  [667,] 0.06090899 0.08124197 0.8578490
##  [668,] 0.06090899 0.08124197 0.8578490
##  [669,] 0.06090899 0.08124197 0.8578490
##  [670,] 0.06090899 0.08124197 0.8578490
##  [671,] 0.06090899 0.08124197 0.8578490
##  [672,] 0.06090899 0.08124197 0.8578490
##  [673,] 0.06090899 0.08124197 0.8578490
##  [674,] 0.06090899 0.08124197 0.8578490
##  [675,] 0.06090899 0.08124197 0.8578490
##  [676,] 0.06090899 0.08124197 0.8578490
##  [677,] 0.06090899 0.08124197 0.8578490
##  [678,] 0.06090899 0.08124197 0.8578490
##  [679,] 0.06090899 0.08124197 0.8578490
##  [680,] 0.06090899 0.08124197 0.8578490
##  [681,] 0.06090899 0.08124197 0.8578490
##  [682,] 0.06090899 0.08124197 0.8578490
##  [683,] 0.06090899 0.08124197 0.8578490
##  [684,] 0.06090899 0.08124197 0.8578490
##  [685,] 0.06090899 0.08124197 0.8578490
##  [686,] 0.06090899 0.08124197 0.8578490
##  [687,] 0.06090899 0.08124197 0.8578490
##  [688,] 0.06090899 0.08124197 0.8578490
##  [689,] 0.06090899 0.08124197 0.8578490
##  [690,] 0.06090899 0.08124197 0.8578490
##  [691,] 0.06090899 0.08124197 0.8578490
##  [692,] 0.06090899 0.08124197 0.8578490
##  [693,] 0.06090899 0.08124197 0.8578490
##  [694,] 0.06090899 0.08124197 0.8578490
##  [695,] 0.06090899 0.08124197 0.8578490
##  [696,] 0.06090899 0.08124197 0.8578490
##  [697,] 0.06090899 0.08124197 0.8578490
##  [698,] 0.06090899 0.08124197 0.8578490
##  [699,] 0.06090899 0.08124197 0.8578490
##  [700,] 0.06090899 0.08124197 0.8578490
##  [701,] 0.06090899 0.08124197 0.8578490
##  [702,] 0.06090899 0.08124197 0.8578490
##  [703,] 0.06090899 0.08124197 0.8578490
##  [704,] 0.06090899 0.08124197 0.8578490
##  [705,] 0.06090899 0.08124197 0.8578490
##  [706,] 0.06090899 0.08124197 0.8578490
##  [707,] 0.06090899 0.08124197 0.8578490
##  [708,] 0.06090899 0.08124197 0.8578490
##  [709,] 0.06090899 0.08124197 0.8578490
##  [710,] 0.06090899 0.08124197 0.8578490
##  [711,] 0.06090899 0.08124197 0.8578490
##  [712,] 0.06090899 0.08124197 0.8578490
##  [713,] 0.06090899 0.08124197 0.8578490
##  [714,] 0.06090899 0.08124197 0.8578490
##  [715,] 0.06090899 0.08124197 0.8578490
##  [716,] 0.06090899 0.08124197 0.8578490
##  [717,] 0.06090899 0.08124197 0.8578490
##  [718,] 0.06090899 0.08124197 0.8578490
##  [719,] 0.06090899 0.08124197 0.8578490
##  [720,] 0.06090899 0.08124197 0.8578490
##  [721,] 0.06090899 0.08124197 0.8578490
##  [722,] 0.06090899 0.08124197 0.8578490
##  [723,] 0.06090899 0.08124197 0.8578490
##  [724,] 0.06090899 0.08124197 0.8578490
##  [725,] 0.06090899 0.08124197 0.8578490
##  [726,] 0.06090899 0.08124197 0.8578490
##  [727,] 0.06090899 0.08124197 0.8578490
##  [728,] 0.06090899 0.08124197 0.8578490
##  [729,] 0.06090899 0.08124197 0.8578490
##  [730,] 0.06090899 0.08124197 0.8578490
##  [731,] 0.06090899 0.08124197 0.8578490
##  [732,] 0.06090899 0.08124197 0.8578490
##  [733,] 0.06090899 0.08124197 0.8578490
##  [734,] 0.06090899 0.08124197 0.8578490
##  [735,] 0.06090899 0.08124197 0.8578490
##  [736,] 0.06090899 0.08124197 0.8578490
##  [737,] 0.06090899 0.08124197 0.8578490
##  [738,] 0.06090899 0.08124197 0.8578490
##  [739,] 0.06090899 0.08124197 0.8578490
##  [740,] 0.06090899 0.08124197 0.8578490
##  [741,] 0.06090899 0.08124197 0.8578490
##  [742,] 0.06090899 0.08124197 0.8578490
##  [743,] 0.06090899 0.08124197 0.8578490
##  [744,] 0.06090899 0.08124197 0.8578490
##  [745,] 0.06090899 0.08124197 0.8578490
##  [746,] 0.06090899 0.08124197 0.8578490
##  [747,] 0.06090899 0.08124197 0.8578490
##  [748,] 0.06090899 0.08124197 0.8578490
##  [749,] 0.06090899 0.08124197 0.8578490
##  [750,] 0.06090899 0.08124197 0.8578490
##  [751,] 0.06090899 0.08124197 0.8578490
##  [752,] 0.06090899 0.08124197 0.8578490
##  [753,] 0.06090899 0.08124197 0.8578490
##  [754,] 0.06090899 0.08124197 0.8578490
##  [755,] 0.06090899 0.08124197 0.8578490
##  [756,] 0.06090899 0.08124197 0.8578490
##  [757,] 0.06090899 0.08124197 0.8578490
##  [758,] 0.06090899 0.08124197 0.8578490
##  [759,] 0.06090899 0.08124197 0.8578490
##  [760,] 0.06090899 0.08124197 0.8578490
##  [761,] 0.06090899 0.08124197 0.8578490
##  [762,] 0.06090899 0.08124197 0.8578490
##  [763,] 0.06090899 0.08124197 0.8578490
##  [764,] 0.06090899 0.08124197 0.8578490
##  [765,] 0.06090899 0.08124197 0.8578490
##  [766,] 0.06090899 0.08124197 0.8578490
##  [767,] 0.06090899 0.08124197 0.8578490
##  [768,] 0.06090899 0.08124197 0.8578490
##  [769,] 0.06090899 0.08124197 0.8578490
##  [770,] 0.06090899 0.08124197 0.8578490
##  [771,] 0.06090899 0.08124197 0.8578490
##  [772,] 0.06090899 0.08124197 0.8578490
##  [773,] 0.06090899 0.08124197 0.8578490
##  [774,] 0.06090899 0.08124197 0.8578490
##  [775,] 0.06090899 0.08124197 0.8578490
##  [776,] 0.06090899 0.08124197 0.8578490
##  [777,] 0.06090899 0.08124197 0.8578490
##  [778,] 0.06090899 0.08124197 0.8578490
##  [779,] 0.06090899 0.08124197 0.8578490
##  [780,] 0.06090899 0.08124197 0.8578490
##  [781,] 0.06090899 0.08124197 0.8578490
##  [782,] 0.06090899 0.08124197 0.8578490
##  [783,] 0.06090899 0.08124197 0.8578490
##  [784,] 0.06090899 0.08124197 0.8578490
##  [785,] 0.06090899 0.08124197 0.8578490
##  [786,] 0.06090899 0.08124197 0.8578490
##  [787,] 0.06090899 0.08124197 0.8578490
##  [788,] 0.06090899 0.08124197 0.8578490
##  [789,] 0.06090899 0.08124197 0.8578490
##  [790,] 0.06090899 0.08124197 0.8578490
##  [791,] 0.06090899 0.08124197 0.8578490
##  [792,] 0.06090899 0.08124197 0.8578490
##  [793,] 0.06090899 0.08124197 0.8578490
##  [794,] 0.06090899 0.08124197 0.8578490
##  [795,] 0.06090899 0.08124197 0.8578490
##  [796,] 0.06090899 0.08124197 0.8578490
##  [797,] 0.06090899 0.08124197 0.8578490
##  [798,] 0.06090899 0.08124197 0.8578490
##  [799,] 0.06090899 0.08124197 0.8578490
##  [800,] 0.06090899 0.08124197 0.8578490
##  [801,] 0.06090899 0.08124197 0.8578490
##  [802,] 0.06090899 0.08124197 0.8578490
##  [803,] 0.06090899 0.08124197 0.8578490
##  [804,] 0.06090899 0.08124197 0.8578490
##  [805,] 0.06090899 0.08124197 0.8578490
##  [806,] 0.06090899 0.08124197 0.8578490
##  [807,] 0.06090899 0.08124197 0.8578490
##  [808,] 0.06090899 0.08124197 0.8578490
##  [809,] 0.06090899 0.08124197 0.8578490
##  [810,] 0.06090899 0.08124197 0.8578490
##  [811,] 0.06090899 0.08124197 0.8578490
##  [812,] 0.06090899 0.08124197 0.8578490
##  [813,] 0.06090899 0.08124197 0.8578490
##  [814,] 0.06090899 0.08124197 0.8578490
##  [815,] 0.06090899 0.08124197 0.8578490
##  [816,] 0.06090899 0.08124197 0.8578490
##  [817,] 0.06090899 0.08124197 0.8578490
##  [818,] 0.06090899 0.08124197 0.8578490
##  [819,] 0.06090899 0.08124197 0.8578490
##  [820,] 0.06090899 0.08124197 0.8578490
##  [821,] 0.06090899 0.08124197 0.8578490
##  [822,] 0.06090899 0.08124197 0.8578490
##  [823,] 0.06090899 0.08124197 0.8578490
##  [824,] 0.06090899 0.08124197 0.8578490
##  [825,] 0.06090899 0.08124197 0.8578490
##  [826,] 0.06090899 0.08124197 0.8578490
##  [827,] 0.06090899 0.08124197 0.8578490
##  [828,] 0.06090899 0.08124197 0.8578490
##  [829,] 0.06090899 0.08124197 0.8578490
##  [830,] 0.06090899 0.08124197 0.8578490
##  [831,] 0.06090899 0.08124197 0.8578490
##  [832,] 0.06090899 0.08124197 0.8578490
##  [833,] 0.06090899 0.08124197 0.8578490
##  [834,] 0.06090899 0.08124197 0.8578490
##  [835,] 0.06090899 0.08124197 0.8578490
##  [836,] 0.06090899 0.08124197 0.8578490
##  [837,] 0.06090899 0.08124197 0.8578490
##  [838,] 0.06902629 0.08931849 0.8416552
##  [839,] 0.06902629 0.08931849 0.8416552
##  [840,] 0.06902629 0.08931849 0.8416552
##  [841,] 0.06902629 0.08931849 0.8416552
##  [842,] 0.06902629 0.08931849 0.8416552
##  [843,] 0.06902629 0.08931849 0.8416552
##  [844,] 0.06902629 0.08931849 0.8416552
##  [845,] 0.06902629 0.08931849 0.8416552
##  [846,] 0.06902629 0.08931849 0.8416552
##  [847,] 0.06902629 0.08931849 0.8416552
##  [848,] 0.06902629 0.08931849 0.8416552
##  [849,] 0.06902629 0.08931849 0.8416552
##  [850,] 0.06902629 0.08931849 0.8416552
##  [851,] 0.06902629 0.08931849 0.8416552
##  [852,] 0.06902629 0.08931849 0.8416552
##  [853,] 0.06902629 0.08931849 0.8416552
##  [854,] 0.06902629 0.08931849 0.8416552
##  [855,] 0.06902629 0.08931849 0.8416552
##  [856,] 0.06902629 0.08931849 0.8416552
##  [857,] 0.06902629 0.08931849 0.8416552
##  [858,] 0.06902629 0.08931849 0.8416552
##  [859,] 0.06902629 0.08931849 0.8416552
##  [860,] 0.06902629 0.08931849 0.8416552
##  [861,] 0.06902629 0.08931849 0.8416552
##  [862,] 0.06902629 0.08931849 0.8416552
##  [863,] 0.06902629 0.08931849 0.8416552
##  [864,] 0.06902629 0.08931849 0.8416552
##  [865,] 0.06902629 0.08931849 0.8416552
##  [866,] 0.06902629 0.08931849 0.8416552
##  [867,] 0.06902629 0.08931849 0.8416552
##  [868,] 0.06902629 0.08931849 0.8416552
##  [869,] 0.07142407 0.09444182 0.8341341
##  [870,] 0.07142407 0.09444182 0.8341341
##  [871,] 0.07142407 0.09444182 0.8341341
##  [872,] 0.07142407 0.09444182 0.8341341
##  [873,] 0.07142407 0.09444182 0.8341341
##  [874,] 0.07142407 0.09444182 0.8341341
##  [875,] 0.07142407 0.09444182 0.8341341
##  [876,] 0.07142407 0.09444182 0.8341341
##  [877,] 0.07142407 0.09444182 0.8341341
##  [878,] 0.07142407 0.09444182 0.8341341
##  [879,] 0.07142407 0.09444182 0.8341341
##  [880,] 0.07142407 0.09444182 0.8341341
##  [881,] 0.07142407 0.09444182 0.8341341
##  [882,] 0.07142407 0.09444182 0.8341341
##  [883,] 0.07142407 0.09444182 0.8341341
##  [884,] 0.07142407 0.09444182 0.8341341
##  [885,] 0.07142407 0.09444182 0.8341341
##  [886,] 0.07142407 0.09444182 0.8341341
##  [887,] 0.07142407 0.09444182 0.8341341
##  [888,] 0.07142407 0.09444182 0.8341341
##  [889,] 0.07142407 0.09444182 0.8341341
##  [890,] 0.07142407 0.09444182 0.8341341
##  [891,] 0.07142407 0.09444182 0.8341341
##  [892,] 0.07142407 0.09444182 0.8341341
##  [893,] 0.07142407 0.09444182 0.8341341
##  [894,] 0.07142407 0.09444182 0.8341341
##  [895,] 0.07142407 0.09444182 0.8341341
##  [896,] 0.07142407 0.09444182 0.8341341
##  [897,] 0.07142407 0.09444182 0.8341341
##  [898,] 0.07142407 0.09444182 0.8341341
##  [899,] 0.07427120 0.09058374 0.8351451
##  [900,] 0.07427120 0.09058374 0.8351451
##  [901,] 0.07427120 0.09058374 0.8351451
##  [902,] 0.07427120 0.09058374 0.8351451
##  [903,] 0.07427120 0.09058374 0.8351451
##  [904,] 0.07427120 0.09058374 0.8351451
##  [905,] 0.07427120 0.09058374 0.8351451
##  [906,] 0.07427120 0.09058374 0.8351451
##  [907,] 0.07427120 0.09058374 0.8351451
##  [908,] 0.07427120 0.09058374 0.8351451
##  [909,] 0.07427120 0.09058374 0.8351451
##  [910,] 0.07427120 0.09058374 0.8351451
##  [911,] 0.07427120 0.09058374 0.8351451
##  [912,] 0.07427120 0.09058374 0.8351451
##  [913,] 0.07427120 0.09058374 0.8351451
##  [914,] 0.07427120 0.09058374 0.8351451
##  [915,] 0.07427120 0.09058374 0.8351451
##  [916,] 0.07427120 0.09058374 0.8351451
##  [917,] 0.07427120 0.09058374 0.8351451
##  [918,] 0.07427120 0.09058374 0.8351451
##  [919,] 0.07427120 0.09058374 0.8351451
##  [920,] 0.07427120 0.09058374 0.8351451
##  [921,] 0.07427120 0.09058374 0.8351451
##  [922,] 0.07427120 0.09058374 0.8351451
##  [923,] 0.07427120 0.09058374 0.8351451
##  [924,] 0.07427120 0.09058374 0.8351451
##  [925,] 0.07427120 0.09058374 0.8351451
##  [926,] 0.07427120 0.09058374 0.8351451
##  [927,] 0.07427120 0.09058374 0.8351451
##  [928,] 0.07427120 0.09058374 0.8351451
##  [929,] 0.07427120 0.09058374 0.8351451
##  [930,] 0.07427120 0.09058374 0.8351451
##  [931,] 0.07427120 0.09058374 0.8351451
##  [932,] 0.07427120 0.09058374 0.8351451
##  [933,] 0.07427120 0.09058374 0.8351451
##  [934,] 0.07427120 0.09058374 0.8351451
##  [935,] 0.07427120 0.09058374 0.8351451
##  [936,] 0.07427120 0.09058374 0.8351451
##  [937,] 0.07427120 0.09058374 0.8351451
##  [938,] 0.07427120 0.09058374 0.8351451
##  [939,] 0.07427120 0.09058374 0.8351451
##  [940,] 0.07427120 0.09058374 0.8351451
##  [941,] 0.07427120 0.09058374 0.8351451
##  [942,] 0.07427120 0.09058374 0.8351451
##  [943,] 0.07427120 0.09058374 0.8351451
##  [944,] 0.07427120 0.09058374 0.8351451
##  [945,] 0.07427120 0.09058374 0.8351451
##  [946,] 0.07427120 0.09058374 0.8351451
##  [947,] 0.07427120 0.09058374 0.8351451
##  [948,] 0.07427120 0.09058374 0.8351451
##  [949,] 0.07427120 0.09058374 0.8351451
##  [950,] 0.07427120 0.09058374 0.8351451
##  [951,] 0.07427120 0.09058374 0.8351451
##  [952,] 0.07427120 0.09058374 0.8351451
##  [953,] 0.07427120 0.09058374 0.8351451
##  [954,] 0.07427120 0.09058374 0.8351451
##  [955,] 0.07427120 0.09058374 0.8351451
##  [956,] 0.07427120 0.09058374 0.8351451
##  [957,] 0.07427120 0.09058374 0.8351451
##  [958,] 0.07427120 0.09058374 0.8351451
##  [959,] 0.07427120 0.09058374 0.8351451
##  [960,] 0.07427120 0.09058374 0.8351451
##  [961,] 0.07427120 0.09058374 0.8351451
##  [962,] 0.07427120 0.09058374 0.8351451
##  [963,] 0.07427120 0.09058374 0.8351451
##  [964,] 0.07427120 0.09058374 0.8351451
##  [965,] 0.07427120 0.09058374 0.8351451
##  [966,] 0.07427120 0.09058374 0.8351451
##  [967,] 0.07427120 0.09058374 0.8351451
##  [968,] 0.07427120 0.09058374 0.8351451
##  [969,] 0.07427120 0.09058374 0.8351451
##  [970,] 0.07427120 0.09058374 0.8351451
##  [971,] 0.07427120 0.09058374 0.8351451
##  [972,] 0.07427120 0.09058374 0.8351451
##  [973,] 0.07427120 0.09058374 0.8351451
##  [974,] 0.07427120 0.09058374 0.8351451
##  [975,] 0.07427120 0.09058374 0.8351451
##  [976,] 0.07427120 0.09058374 0.8351451
##  [977,] 0.07427120 0.09058374 0.8351451
##  [978,] 0.07427120 0.09058374 0.8351451
##  [979,] 0.07427120 0.09058374 0.8351451
##  [980,] 0.07427120 0.09058374 0.8351451
##  [981,] 0.07427120 0.09058374 0.8351451
##  [982,] 0.07427120 0.09058374 0.8351451
##  [983,] 0.07427120 0.09058374 0.8351451
##  [984,] 0.07427120 0.09058374 0.8351451
##  [985,] 0.07427120 0.09058374 0.8351451
##  [986,] 0.07427120 0.09058374 0.8351451
##  [987,] 0.07427120 0.09058374 0.8351451
##  [988,] 0.07427120 0.09058374 0.8351451
##  [989,] 0.07427120 0.09058374 0.8351451
##  [990,] 0.07427120 0.09058374 0.8351451
##  [991,] 0.07427120 0.09058374 0.8351451
##  [992,] 0.07427120 0.09058374 0.8351451
##  [993,] 0.07427120 0.09058374 0.8351451
##  [994,] 0.07427120 0.09058374 0.8351451
##  [995,] 0.07427120 0.09058374 0.8351451
##  [996,] 0.07427120 0.09058374 0.8351451
##  [997,] 0.07427120 0.09058374 0.8351451
##  [998,] 0.07427120 0.09058374 0.8351451
##  [999,] 0.07427120 0.09058374 0.8351451
## [1000,] 0.07427120 0.09058374 0.8351451
## [1001,] 0.07427120 0.09058374 0.8351451
## [1002,] 0.07427120 0.09058374 0.8351451
## [1003,] 0.07427120 0.09058374 0.8351451
## [1004,] 0.07427120 0.09058374 0.8351451
## [1005,] 0.07427120 0.09058374 0.8351451
## [1006,] 0.07427120 0.09058374 0.8351451
## [1007,] 0.07427120 0.09058374 0.8351451
## [1008,] 0.07427120 0.09058374 0.8351451
## [1009,] 0.07427120 0.09058374 0.8351451
## [1010,] 0.07427120 0.09058374 0.8351451
## [1011,] 0.07427120 0.09058374 0.8351451
## [1012,] 0.07427120 0.09058374 0.8351451
## [1013,] 0.07427120 0.09058374 0.8351451
## [1014,] 0.07427120 0.09058374 0.8351451
##      Expert_1  Expert_2  Expert_3
## T1  0.3333333 0.3333333 0.3333333
## T2  0.3333333 0.3333333 0.3333333
## T3  0.3333333 0.3333333 0.3333333
## T4  0.3333333 0.3333333 0.3333333
## T5  0.3333333 0.3333333 0.3333333
## T6  0.3333333 0.3333333 0.3333333
## T7  0.3333333 0.3333333 0.3333333
## T8  0.3333333 0.3333333 0.3333333
## T9  0.3333333 0.3333333 0.3333333
## T10 0.3333333 0.3333333 0.3333333
myT <- length(returns)
T0  <- 1250   
K   <- length(gammaGrid)
print(K)
## [1] 3
date <- index(returns)[T0 : (myT - 1)]


##############################################

err <- results$err
errSeqOC <- results$errSeqOC
print(err)
##    [1] 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0
##   [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##   [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [112] 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [149] 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0
##  [186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
##  [223] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
##  [260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [297] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [334] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [371] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
##  [408] 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
##  [445] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
##  [482] 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [519] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [556] 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0
##  [593] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [630] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0
##  [667] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [704] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
##  [741] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [778] 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
##  [815] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 1
##  [852] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [889] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
##  [926] 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [963] 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## [1000] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
mean(err)
## [1] 0.06213018
for (k in 1:K) {
  print(mean(errSeqOC[k,-1]))
}
## [1] 0.09772952
## [1] 0.100691
## [1] 0.1026654
par(mfrow = c(1, 1)) 


colors <- rainbow(K)


plot(date[-1], cummean(err[-1]), type = "l", col = "black",
     ylim = c(0, 0.5),
     main = "Cumulative Errors",
     xlab = "Data", ylab = expression(alpha))


abline(h = 0.1, col = "red", lwd = 0.3, lty = 2)


for (k in 1:K) {
  lines(date[-1], cummean(errSeqOC[k,-1]), col = colors[k], lwd = 1.5)
}

#Legend
legend("topright", legend = c("Aggregated", paste("Expert", 1:K)),
       col = c("black", colors), lty = 1, lwd = 1.5, cex = 0.8)

In the first plot, the cumulative errors of the three experts stabilize over the long term around the target value of 0.1, showing convergence in their performance. However, the aggregated predictor consistently stays at lower levels, suggesting it is able to reduce errors more effectively than the individual experts, likely due to the combination of information from all models.

# LCR across gammas
lcr_best_gamma_05 <- compute_local_coverage(errSeqOC_EStM, window = 250)  # Gamma = 0.05
lcr_best_gamma_03 <- compute_local_coverage(errSeqOC_St_003, window = 250)  # Gamma = 0.03
lcr_best_gamma_07 <- compute_local_coverage(errSeqOC_EStM_007, window = 250)  # Gamma = 0.07
lcr_aggregated <- compute_local_coverage(err, window = 250)
valid_idx <- !is.na(lcr_aggregated)
date_valid <- date[valid_idx]
lcr_aggregated_valid <- lcr_aggregated[valid_idx]
lcr_best_gamma_05_valid <- lcr_best_gamma_05[valid_idx]
lcr_best_gamma_03_valid <- lcr_best_gamma_03[valid_idx]
lcr_best_gamma_07_valid <- lcr_best_gamma_07[valid_idx]
colors_best <- c("blue", "green", "purple")

# LCR
plot(date_valid, lcr_aggregated_valid, type = "l", col = "black",
     ylim = c(0.8, 1), main = "Comparison Local Coverage Rate (LCR)",
     xlab = "Date", ylab = "LCR")


abline(h = 0.9, col = "red", lwd = 0.3, lty = 2)


lines(date_valid, lcr_best_gamma_05_valid, col = colors_best[1], lwd = 1.5)
lines(date_valid, lcr_best_gamma_03_valid, col = colors_best[2], lwd = 1.5)
lines(date_valid, lcr_best_gamma_07_valid, col = colors_best[3], lwd = 1.5)

# Legend
legend("bottomright", legend = c("Aggregated", "Gamma = 0.05", "Gamma = 0.03", "Gamma = 0.07"), 
       col = c("black", colors_best), lty = 1, lwd = 1.5, cex = 0.8)

In the second plot, the Local Coverage Rates (LCR) of the three models with different gamma values fluctuate similarly around the target value of 0.9. The aggregated predictor, however, consistently stays at higher levels, often above 0.95, converging towards the target only at the end. This could indicate that the aggregated model generates wider intervals compared to the individual models, thereby including a larger range of values but deviating from the target. This is neither necessarily positive nor negative: it could suggest greater reliability or, conversely, a tendency to overcover, sacrificing accuracy.

Final Choice

Among the models analyzed, we would choose the Adaptive Momentum eGARCH with gamma = 0.07 as it achieves a local coverage rate (LCR) closer to the target value. The aggregated model, on the other hand, shows a significantly higher LCR, reaching levels around 0.95, which is well above the desired target. This makes the Adaptive Momentum eGARCH with gamma = 0.07 a more balanced and reliable option for practical use, ensuring better alignment with the target coverage rate.